A Search for Cosmic Microwave Background Anisotropies on 
Arcminute Scales with Bolocam 
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ABSTRACT 



' We have surveyed two science fields totaling one square degree with Bolocam 

at 2.1 mm to search for secondary CMB anisotropies caused by the Sunyaev- 



^ ■ Zel'dovich effect (SZE). The fields are in the Lynx and Subaru/XMM SDSl 

fields. Our survey is sensitive to angular scales with an effective angular multi- 
^ • pole of ieff = 5700 with FWHM^ = 2800 and has an angular resolution of 60 

^ ■ arcseconds FWHM. Our data provide no evidence for anisotropy. We are able 

■ to constrain the level of total astronomical anisotropy, modeled as a fiat band 
^ ' power in with frequentist 68%, 90%, and 95% CL upper limits of 590, 760, 

and 830 /ii^^^j^. We statistically subtract the known contribution from pri- 
CN ■ mary CMB anisotropy, including cosmic variance, to obtain constraints on the 

SZE anisotropy contribution. Now including flux calibration uncertainty, our 

■ frequentist 68%, 90% and 95% CL upper hmits on a fiat band power in Ce are 
Cf) • 690, 960, an d 1000 ^K%Mpt. Wh e n we instead employ the analytic spectrum 



lO ' suggested by iKomatsu and Seljakl (120021 ) . and account for the non-Gaussianity 



^ ■ of the SZE anisotropy signal, we obtain upper limits on the average amplitude of 

O ■ their spectrum weighted by our transfer function of 790, 1060, and 1080 /ii^^^^j^. 

We obtain a 90% CL upper limit on ag, which normalizes the power spectrum 
^ ■ of density fluctuations, of 1.57. These are the first constraints on anisotropy and 

■ (Tg from survey data at these angular scales at frequencies near 150 GHz. 
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Subject headings: cosmology: observation — cosmic microwave background — 
methods: data analysis — large-scale structure of the universe — cosmological 
parameters 



Introduction 



1.1. Background 

The SZeEI is the inverse Compton scattering of CMB ph otons with a distribution of ho t 
electrons, causing a net increase in the energy of the photons fISunyaev and Zel'dovichlll972l ). 
Since the background CMB is redshifted along with the SZE-induced distortion, the relative 
amplitude of the distortion, ATcmb/Tcmb, is independent of redshift. The distortion caused 
by the SZE is proportional to the Comptonization parameter y, which is a measure of the 
integral of the electron thermal energy density along the line of sight and is given by 



y 



dl TleksTe, 



where o"r is the Thomson cross section, rrie is the electron mass, c is the speed of light, is 
Boltzmann's constant, Te is the temperature of the electrons, and rig is the number density 
of electrons. Since the scattering process conserves photon number, the thermal spectrum 
of the CMB is distorted by the SZE; there is a negative temperature shift at low frequency 
and a positive temperature shift at high frequency. The cross-over point where there is no 
distortion of the CMB occurs at approximately 218 GHz. The temperature shift caused by 
the SZE, ATcMB, is 

ATcMB 



-CMB 



f{x)y, 



where 



fix) 



X - 4, 



— 1 

and X = hu/ksTcMB, h is Planck's constant, u is the frequency, and Tcmb = 2.73 K is the 
temperature of the CM B. For reference, excell ent reviews of the SZE and its relevance to 
cosmology are given by iBirkinshawi (1l999l ) and ICarlstrom et al.l (l2002l ) . 



^Throughout this paper SZE refers to the thermal SZE. 
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1.2. Untargeted SZE Surveys 



To date, there have been no detections of previously unknown clusters using the SZE. 
However, unresolved objects in SZE surveys will produce anisotropies in the CMB that are 
expected to dominate the CMB power spectrum at small angular scales corresponding to 
angular multipoles above I ~ 2500. The overall normalization of these SZE-induced CMB 
anisotropies is extremel y sensitive to cts and can be used to constrain the value of this 
cosmological parameter (IKomatsu and Seljakl |2002| ) . Several experiments have conducted 
SZE surveys that have produced tentative detections of the SZE-induced anisotropies in 
the CMB. At 30 GHz CBI has me asured an excess C MB power between C. = 2000 and 
i = 3500 at a significance of 3. la jMason et al.ll2003h . Also at 30 GHz, BIMA/OVRO 
has measured a CM B anisotropy of 220]'^}2o I^^cmb angular multipole oi i = 5237 

dPawson et aPbood ). AC BAR, at 150 GHz and 2000 < i < 3000, has measured an excess 
power of 34 ± 20 yuK^^^ flReichardt et allboosh . A joint analysis of the CBI and ACBAR 
excesses sho ws that they are six tim es more likely to be caused by the SZE than primordial 
fluctuations (IReichardt et al.ll2008l ). 



Additionally, these tentative anisotropy detections have be en used to constra in cosmo- 
logical parameters; the CBI data are consist ent with as ^ 1 (IBond et al.l 120051). and the 
BIMA/OVRO data measure ag = 1.03;|];^^ dPawson et aPbooeh . IReichardt et all J2008h 
combine various data sets to place constraints on ag via an excess contribution to anisotropy 
at high i, i > 1950. ACBAR and WMAP3, plus CBI and BIMA/OVRO data at high i 
and lower frequency, combine to indicate af^ = 0.95^*^0 04 when the amplitude of the SZE 
contribution is not slaved to its contribution at low i. Since this result is inconsistent with 
other constraints on g g, including those from the lower i portions of the CMB spectrum, 
Reichardt et al.l (120081 ) also consider a case in which the high-£ SZE contribution is slaved to 



to the contribution at lower i via Ug = erg. In this case the excess power is produced by point 
sources. Fitted to the CMBall data set, which excludes CBI high-£ and BIMA/OVRO data, 
this model results in erg values consistent with other measurements, ~ 0.80 — 0.81±0. 03 — 0.04 
depending on the assumptions and data sets included. No attempt is made to explain the 
CBI and BIMA/OVRO excesses. Overall, the current results suggest two possibilities: there 
are point source contributions to all the high-£ data (ACBAR, CBI, BIMA/OVRO) that have 
not been properly included; or the SZ contribution calculated from theory is underestimated. 

The survey presented here is the first such survey at 150 GHz and at ^ ~ 6000. As we 
shall explain, contributions from primary CMB anisotropies, SZE, radio, and submillimeter 
point sources are all expected to be comparable, each at a level of Ce ~ 50 /iK^j^^^. 
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Observations 



2.1. Instrument Description 

Bolocam is a large format, 144 detector, millimeter- wave camera designed to be op- 
erated at the Caltech Submillimeter Observatory (CSO). For these observations the array 
was comprised of 115 optical and 6 dark detectors. Each detector is housed within its 



own integrating cavity, formed by a frontshort plate and a backshort plate (IGlenn et al. 



2OO2I ) ■ Smooth-walled conical feedhorns separated by 0.7 (f/#)A, a cold (4 K) high-density 
polyethylene (HDPE) lens, and a room-temperature ellipsoidal mirror are used to couple the 
detectors to the CSO optics. Each feedhorn terminates into a cylindrical waveguide, which 
defines the low-frequency cutoff of the system; the final filter in a series of six cold metal- 
mesh filters determines the high-frequency cutoff. The resulting passband is centered at 
143 GHz, and has an effective width of 21 GHz. A cold (4 K) Lyot stop is used to define the 
illumination of the 10.4 m primary mirror, and the resulting far-field beams have FWHMs 
of 60 arcseconds. Bolocam can also observe at 270 GHz, and has been used in this mode for 



(Laurent et al. 


2005; 


Enoch et al. 


2006; 


Young et al. 


2006) 



The detector array has a hexagonal georn etry, and utilizes silicon nitride micromesh 
(spider-web) bolometers ( Mauskopf et al. 1997) which are cooled to 260 mK using a three- 
stage ^He/^He/^He sorption refrigerator flBhatia et al.ll200d . l2002[ ). JFETs located near the 
array and operated at 140 K are used to buffer the high-impedance bolometer signals from 
sources of current noise. In order to avoid the 1// noise from the JFETs, the bolometers are 
biased at 130 Hz and read out using room-temperature loc k in amplifiers . More details of the 



Boloc am in strument can be fo und in iGolwala et al.l (120081 ) , iGlenn et al.l (119981 ) , iGlenn et al. 
(120031 ) ■ and lHaig et aP tooi ). 



2.2. Observing Strategy 

The data described in this paper was collected during a forty night observing run in late 
2003. During the first half of each night we observed a 0.5 deg^ region centered at 02hl8m00s, 
-5d00m00s (J2000), which coincides with the Subaru/XMM Deep Survey (SXDS or SDSl); 
and during the second half of each night we observed a 0.5 deg^ region centered on the 
Lynx field at 08h49ml2s, +44d50m24s (J2000). These fields were selected because they 
have extremely low dust emission and a large amount of optical/X-ray data that could be 
used to follow up any SZE cluster candidates found in the maps. 
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Two nights at the start of the run were used to analyze different scan strategies for 
mapping the science fields. The maps were made by repeatedly raster scanning across the 
field, stepping perpendicular to the scan, then rastering across the field in the opposite 
direction until the entire field has been covered. Our studies showed that the time-stream 
noise is independent of the angle of the raster scan and the turnaround time between scans, 
so we chose to scan parallel to RA or dec and turnaround as quickly as the telescope would 
allow (~ 10 seconds). Additionally, we found that our sensitivity to astronomical signals 
is maximized when we raster scan at a speed of 240 arcseconds/second|§ At this speed it 
takes approximately 12.5 seconds to complete one scan across the field, which means we were 
on-source approximately 56% of the time during an observation. Although scanning at this 
relatively quick speed reduces our time on-source because a larger fraction of time is spent 
on turnarounds between scans, it also puts a larger amount of our signal band above the 
1// atmospheric noise. Given the scan speed and turnaround time mentioned above, along 
with our step size of 162 arcseconds (~ 1/3 of the field of view) a complete map of the field 
was made in approximately eight minutes. 



3. Data Reduction 

3.1. Initial Processing 

After merging the bolometer time-streams recorded by the data acquisition system with 
the pointing information recorded by the telescope, we parse the data into files that contain 
a single observation. Each single observation contains a set of scans that completely map the 
astronomical field or object, and they are typically around ten minutes in length. Parsing the 
data by observation is useful because individual observations are statistically independent, 
have a small enough number of data samples to be easily manageable from an analysis 
standpoint, and provide a convenient division of the data for the sake of bookkeeping. 

Once the initial merging and parsing of the data is complete, we begin the process of 
refining the data. The first step in this process is to remove the effects of the lockin amplifier 
electronics filters and to down-sample the data from 50 Hz to 10 Hz. We down-sample the 
data because essentially no astronomical signal is lost, while a large amount of 60- Hz pickup 
noise is removed. See Figured] for an illustration of the noise spectrum and the shape of the 
beam in frequency space. 



^ Faster speeds were not attempted due to fears that the CSO would not function properly and/or would 
be damaged. 
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3.2. Noise Removal 

There are several forms of correlated noise present in the raw bolometer data which 
contaminate the astronomical signal and therefore must be modeled and removed. First, 
the emission from the atmosphere changes as a function of telescope elevation angle due to 
the changing path length through the atmosphere. The path length through the atmosphere 
relative to the zenith path length is called the airmass, A, and is described by 

A = l/sin(e), 

where e is the elevation angle. For a typical observation the range of elevation angles is 
approximately one degree, which corresponds to a change in airmass between 0.005 and 0.060 
for elevation angles between 75 and 30 degrees. For reference, a change of 0.060 in airmass 
corresponds to a change of approximately 0.5 K of optical loading from the atmosphere, 
or a change in surface brightness of a little less than 1 Kcmb- To remove this elevation- 
dependent signal, we calculate a linear fit of bolometer signal versus airmass. We build up a 
fit using each 12.5-second-long scan within the observation, after removing the mean signal 
level and airmass for the scan. This process yields one set of linear fit coefficients for each 
bolometer for the entire observation, which is used to create a template that is removed from 
the bolometer time-streams. 

Next, we create a template from the bias voltage monitors to account for the small 
amount of noise from the bias electronics. Note that the bias applied to the bolometers 
is monitored through amplifier electronics identical in design to those used to monitor the 
bolometer signals. This template is then correlated and removed from each of the the 
bolometer time-streams. A template is also created from the dark bolometer signals and 
removed from the bolometer time-streams. Note that both the bias template and dark 
bolometer template have an RMS of < 1 toKcmb- 

Finally, and most importantly, we remove a template describing the fiuctuations in 
emission from the atmosphere [i.e., the atmospheric noise). Since atmospheric noise is the 
dominant signal in our data, and the beams from the individual detectors overlap to a 
high degree while passing through the atmosphere, a template for the atmospheric signal 
is created by averaging the signals from all the bolometersjfl Three different algorithms 
are used to construct this template, one for which the atmospheric signal is assumed to be 
constant over the array, one for which the atmospheric signal is allowed to vary linearly with 



To quantify the degree of overlap, note that the far-field distance for Bolocam at the CSO is approxi- 
mately 30 km; at a height of 20 km, well above most of the water vapor in the atmosphere, the ray bundles 
from nearby detectors have only separated by one half-width. 
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bolometer location on the array, and one for which the atmospheric signal is allowed to vary 
quadratically with bolometer location on the array. 

For the most basic case of an average template, the algorithm proceeds as follows. 
Initially, a template is constructed according to 

V'^^" c d- 



where n is the sample number, is the number of bolometers, Cj is the relative responsivity 
of bolometer z, (ij„ is the signal recorded by bolometer i at time sample number n, and T„ is 
the template. The template generally has an RMS between 10 and 100 mKcMB, depending 
on the observing conditions. A separate template is computed for each 12.5-second-long 
scan. After the template is computed, it is correlated with the signal from each bolometer 
to determine the correlation coefficient, with 



Ej = Ns rp2 
7=1 ^ n 



(2) 



Cj is the correlation coefficient of bolometer i and A^^ is the number of samples in the 12.5- 
second-long scanj^ Next, the q in Equation [1] are set equal to the values of q found from 
Equation [2], and a new template is computed. The process is repeated until the values of 
Cj stabilize. We generally iterate until the average fractional change in the CjS is less than 
1 X 10"*^, which usually takes five to ten iterations. If the qs fail to converge after 100 
iterations, then the scan is discarded from the data. For the more advanced planar and 
quadratic algorithms, the process proceeds in the same way except linear and quadratic 
variations with bolometer position are allowed when the template is constructed. These 
algorithms, along with adaptive PCA an d time-lagged aver age template subtraction, are 



described and compared in more detail in ISayers et al.l (120081 ) 



Each of the three different atmospheric noise removal algorithms, average, planar, and 
quadratic template removal, was applied to each observation. Therefore, three different 
atmospheric-noise-cleaned time-streams are generated for each observation. A figure of merit 
is calculated for each of the three files for each observation, based on the noise level of 
the data and the expected astronomical signal shape. Details of the calculation of this 
figure of merit are given in Section [71 For each observation, the file with the best figure of 
merit value will be the one used to create the final map of the data. Weather is the main 
criteria that determines which algorithm will be selected as optimal for a given observation; 



^ The best fit correlation coefiicients cliange from one scan to tlie next, typically by a couple percent. 
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more aggressive algorithms (planar or quadratic) are selected in poor weather conditions 
and more benign algorithms (average or planar) are selected in good weather conditions. 
However, there is some dependence on the profile of the source, and observations of compact 
objects tend to be optimally processed using more aggressive algorithms than observations 
of extended objects. 

4. Calibration 

4.1. Pointing reconstruction 

Pointing reconstruction consists of determining the location of each detector's beam on 
the sky at each instant in time. We compute this location in two steps: 1) we calculate the 
location of each bolometer relative to the center of the array and 2) we then determine the 
absolute coordinates of the center of the array. 

To determine the relative locations of the bolometers, we observed Uranus or Neptune 
for approximately fifteen minutes every other night. These planets are bright enough to 
appear at high signal-to-noise in a map made from a single bolometer, so they can be used to 
determine the position of each detector relative to the array center. Since Bolocam was held 
at a fixed angle in the alt / az coordinate system for the entire observing run, each bolometer 
views the optics in the same way for the entire run and the coordinates on the sky in alt / az 
units remain fixed. Therefore, we combined the data from all the planet observations to 
determine the average position of each beam on the sky. See Figure O The uncertainties 
on these average positions were ~ 1 arcsecond, which is negligible when compared to the 
60 arcsecond FWHM of a Bolocam beam. We found no evidence for a systematic difference 
in the beam positions derived from any single observation to the average beam position found 
from all the observations. This indicates that the optical system was very stable over the 
entire observing run, including a wide range of telescope elevation angles. 

To determine the absolute location of the center of the array, we observed a bright quasar 
with a known position near the science field for approximately ten minutes once every two 
hours. Three different quasars were used for the SDSl field (0106-1-013, 0113-118, and 0336- 
019), and two different quasars were used for the Lynx field (0804-1-499 and 0923-1-392). Each 
source was observed for five minutes while scanning parallel to RA, then for five minutes 
scanning parallel to dec (analogous to how the science fields were observed). We found no 
systematic offset based on scan direction; the maps made while scanning parallel to RA 
produce the same source location as the maps made while scanning parallel to dec. The 
difference in the centroid location for these consecutive observations was then used to deter- 
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mine the measurement uncertainty for the centroided location of each source. As expected, 
the uncertainty in the centroided location of the five sources is inversely proportional to the 
flux of the source. Additionally, we found no evidence that the measurement uncertainty 
degrades or improves as a function of time during the night for our typical observing times 
between 20:00 and 07:00 local time. 

The pointing data were broken up into three distinct subsets corresponding to the 
azimuthal position of the telescope: SDSl was observed between azimuth angles of 90 and 
270 (in the south), while Lynx was observed between azimuth angles of -90 and 90 (in the 
north approaching from the east), and also between azimuth angles of 270 and 360 (in the 
north approaching from the west). Most of the Lynx data were taken between an azimuth 
angle of -90 to 90, so the third subset of data is considerably smaller than the first two (about 
1/5 the size). Note that the slewing limits of the telescope are roughly equal to azimuth 
angles of -90 to 360. There is a correlation between the elevation angle of the telescope and 
the pointing offset for each of these subsets. We attempted to model this correlation with 
several low-order polynomials, but we found that a quadratic fit of pointing offset versus 
elevation was sufficient since higher-order fits did not significantly reduce the scatter of the 
data. We found no correlation between the telescope azimuth angle and the residual offset, 
other than the slight difference between the pointing models determined for the three subsets. 
Therefore, a simple quadratic fit of pointing offsets versus telescope elevation angle served 
as our only pointing model. Plots of this final model can be found in Figure [3l 

For each of the three subsets, we calculated the uncertainty in the pointing model by 
analyzing the residual offset of each centroid location from the model. Some residual scatter 
is expected due to the measurement uncertainty of each centroid, however the scatter we 
find is slightly larger. The difference between the actual scatter and the predicted scatter 
is consistent for all three subsets, and translates to an uncertainty in the pointing model 
of 4.9 arcseconds. This uncertainty is small compared to our beam size and thus made a 
negligible difference in the beam shape used in the final science analysis. 



4.2. Flux Calibration 



Our flux calibration techniq ue, summarized below , has been used previously with Bolo- 
cam to calibrate 1.1 mm data (ILaurent et al.l l2005l ). Since the amount of astronomical 
signal attenuation by the atmosphere is a function of opacity and airmass, the standard 
flux calibration technique for millimeter-wave instruments requires frequent observations of 
calibration sources that are close to the science field. However, we were able to use a more 
advanced technique with Bolocam because we continuously monitor the operating resistance 
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of the bolometers using the carrier amphtude measured by the bolometer voltage at the 
bias frequency. When the atmospheric transmission decreases, the optical loading from the 
atmosphere increases, which lowers the bolometer resistance. Additionally, the bolometer 
responsivity is a monotonically decreasing function of the bolometer resistance. Therefore, 
by fitting the flux calibration as a function of the bolometer operating resistance, we can 
simultaneously account for changes in the atmospheric transmission and bolometer respon- 
sivity. 

Six different flux calibrations were needed for our data set. The base temperature of 
the sub-Kelvin refrigerator was changed on November 4, 2003, and the bias voltage applied 
to the bolometers was changed on November 5, 8 (twice), and 10, 2003. Each of the bias 
changes caused a change in the responsivity of the bolometers, so a different flux calibration 
is needed after each change. Since the first five data sets are relatively short in duration, 
the observing conditions were relatively constant within each set. Therefore, a constant flux 
calibration, rather than a flux calibration that varies as a function of bolometer operating 
resistance, was adequate to describe the data for these five sets. However, a fit of the flux 
calibration as a function of bolometer operating resistance was required for the final data 
set. 

The relative calibration of the detectors was determined from the science field obser- 
vations. Since these observations covered regions of the sky with negligible amounts of 
astronomical flux, the fluctuations in thermal emission from the atmosphere are the dom- 
inant source of the signal recorded by each bolometer. Additionally, this signal should be 
the only one that is correlated among all the bolometers since the beams from all bolome- 
ters overlap to a high degree when passing through the atmosphere. Therefore, this signal 
should be the same in each bolometer, weighted by the responsivity of that bolometer. So, 
by determining how correlated the data from each bolometer is with this common signal, it 
is possible to determine the relative calibration of each bolometer. The uncertainties in the 
relative calibrations determined using this method are less than 1%. 

The absolute flux calibration was determined from observations of Uranus, Neptune, 
0923+392, and NGC2071IR. Since we did not have enough observations of Uranus and 
Neptune to adequately determine the shape of the calibration versus bolometer operating 
resistance, we used 0923+392 and NGC2071IR as secondary calibrators. These two sources 
are known to have minimal variations in emitted flux as a function of time, so they are well 
suited to be used for determining the functional form of the flux c alibration versus bolometer 



operating resistance relationship (jPeng et al.ll2000l : ISandelllll994l ). Note that we did not use 



any of the published fluxes for 0923+392 or NGC2071IR, rather the fluxes were left as 
free parameters and they were used to determine the shape of the calibration curve versus 
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bolometer operating resistance. We used the peak signal and median bolometer operating 
resistance from each observation to determine the fit parameters in the function 



bolo ) 



bolo J 



where Vj is the peak bolometer signal (in nV) recorded for the j*'^ source, Rboio is the bolome- 
ter operating resistance, Fj is equal to the flux of the j^^ source (known for Uranus and Nep- 
tune, left as a free parameter for NGC2071IR and 0923-1-392), and ai and free parameters. 
The p lan et fluxes were deter mined from the temperature spectra given in iGriffin and Orton 
( I1993I ) or lOrton et al.l (119861 ) , along with the planet solid angles calculated from the planet 
flux calculator at the James Clerk Maxwell Telescope websitejf] For reference, the abso- 
lute calibration ranges from approximately 180 nV/Jy up to 280 nV/Jy over the range of 
bolometer operating resistances recorded during our observing run. See Figure |H 

Our flux calibration uncertainty was determined as follows. First, the temperature pro- 
files of Uranus and Neptune were derived by Griffin and Orton using Mars as an absolute cal- 
ibrator (IGriffin and Ortonlll993l ). To determine the surface brightness of Mars at millimeter 
wavelengths. Griffin and Orton use d the raodel developed by Wright based on observations 
made at far-infrared wavelengths (IWrightl 119761) , alon g with the logarithmic interpolation 
to longer wavelengths described by Griffin et al.l (19861) . The estimated uncertainty on this 
interpolated model is approximately 5% ( Wright 19761 )1^ Second, the uncertainties on the 
temperat ure profiles of Uranus a nd Neptune are estimated to be less than 1.5% relative 
to Mars JOriffin and Ortonlll993h FI Additionally, the observations of Uranus and Neptune 
were taken with a precipitable water vapor of 1.5 ± 0.5 mm, which results in a calibration 
uncertainty of ~ 1.4%. Finally, the error inferred by the scatter of our measurements results 
in calibration uncertainties between 0.6% and 3.0% for each of the data sets. The end result 
is an overall flux calibration uncertainty of approximately 5.5%, limited by the uncertainty 
in the temperature of Mars. 



^ [http : / / www .jach.hawaii.edu/jacbin/ planetflux . pi 

^ There is also a brightness model based on a physical model of the dielectric properties of the Martian 
surface that was developed by Rudy (|Rudvlll987t iRudv et al.lll987l) . This model was constrained by mea- 



surements at centimeter wavelengths, and also needs to be extrapolated to millimeter wavelengths. Grifhn 
and Orton, along with Goldin, et al., compared the results of these two models at millimeter wavelengths , 
and foimd that they agree within their estimated uncertainties ( Griffin and Orto"nlll993t Goldin et al. Il997h . 
Based on the comparison of these two models, Griffin and Orton conclude the the uncertainty in the Martian 
brightness based on the Wright model is 5%. 

^ Grifhn and Orton find that the uncertainty is 1.7 K for both their Uranus and Neptune models. Since 
the temperature of these planets in our band is approximately 115 K, this translates to an uncertainty of 
~ 1.5%. 



4.3. Beam Calibration 



Since the astronomical signals in our maps are inherently smoothed based on the profile 
of the Bolocam beams, it is important to understand their shapes. Additionally, our flux 
calibration is based on observations of point sources, so our maps have units of flux den- 
sity. However, since the CMB or SZE signal we are looking for is a surface brightness or 
temperature, we need to know the area of our beam in solid angle to convert our maps to 
surface brightness units. Therefore, any error in our determination of the beam area will 
show up as a surface brightness or temperature calibration error. To determine the profile of 
our beam, we used the observations of Uranus and Neptune. These planets are well suited 
for measuring our 60 arcsecond FWHM beams; they have semi-diameters of ~ 1 arcsecond, 
which means they are essentially point-like and thus will appear in our maps with shapes 
given by our beam profile. 

Based on simulations, we expected all of the beams to have a similar profile. However, 
we first calculated the beam for each bolometer separately to validate this expectation. There 
was not enough data from a single planet observation to make a high signal-to-noise mea- 
surement of the beam for an individual bolometer, so we averaged the data for groups of four 
bolometers that are close to each other on the focal plane. Nearby bolometers have beams 
with similar paths through the optics, so they should also have similar profiles. Each bolome- 
ter was grouped into four distinct sets, each of which contained four nearby bolometers, and 
the average profile from these four sets was determined. The measurement uncertainty on 
these profiles can be quantified by the standard deviation of the peak-normalized areas of 
the beam profiles, which was approximately 3.1%. See Figure [51 Within our measurement 
uncertainty, all of the individual bolometer beam profiles were consistent, so a single beam 
profile can be used to describe every bolometer. To measure this single beam profile, we 
averaged the data from all of the planet observations for all of the bolometers. The peak- 
normalized area of this profile is 3970 arcseconds^, which is the area of a Gaussian beam 
with a FWHM of 59.2 arcseconds. However, the beam profile is not exactly Gaussian, and 
the measured profile was used for all of our analysis. Since we cannot rule out systematic 
variations in the beam area from one bolometer to the next at the level of our single bolome- 
ter measurement uncertainty, we have conservatively estimated the uncertainty in this beam 
area measurement to be 3.1%. 
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5. Map Making 

5.1. Least Squares Map Making Theory 

The astronomical signals we seek can be thought of as two-dimensional objects, which 
can be represented by a map with finite pixelization. For simplicity, this two-dimensional 
map can be thought of as a vector, m. This map is stored in the bolometer time-streams, d, 
according to 

d = prh + n, (3) 

where p is a matrix containing the pointing information and n is noise. Note that we 
represent matrices with a bold symbol, and vectors with an arrow. Since fh is what we are 
fundamentally interested in obtaining, we need to find a solution to Equation [3] that yields 
the optimum unbiased estimate of rh given d. There are several methods that can be used to 



estimate m, includ ing the commonly used least squares method described below (ITegmark 



19971 : IWrightl 119961 ). 

Solving the least squares problem for Equation [3] requires minimizing 

= (d — prh)^w{d — pm), (4) 

where w is the inverse of the time-stream noise covariance matrix, (nvF^ ^ . The estimator 
for m derived from Equation H] is 

fh! = cp^wd, (5) 

where c = (p"^wp)~^ is the map-space noise covariance matrix. If the time-stream noise, n, 
has a white spectrum, then the various terms in Equation [S] are easy to understand because 
w and c are both diagonal, w is the inverse of the time-stream noise variance, and applies the 
appropriate weight to each sample in the time-stream, p^ then bins the data time-stream 
into a map, and c corrects for the fact that p"^ sums all of the data in a single map bin 
instead of averaging it. The general idea is the same for non- white time-stream noise, but 
w will mix time samples and c will mix map pixels. 

If the time-stream noise is stationary, then the time-stream noise covariance matrix can 
be diagonalized by applying the Fourier transform operator, F. In this case, any element of 
the inverse time-stream noise covariance matrix can be described by 

w(ti,t2) = (n(ii)n{hf) ' = w(At), 



where ti and ^2 are any two time samples separated by At. The corresponding elements of 
the Fourier transform of the inverse covariance matrix, W = FwF~^, can be written as 



W(/i,/2)=W(/i)5/, 



/2' 
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where 5/i,/2 represents a Kronecker delta and / is frequency in The diagonal elements 
of W are equal to 1/(PSD*A/), where PSD is the noise power spectral density and A/ is 
the frequency resolution of the time-stream. The Kronecker delta ensures that all of the 
off-diagonal elements are equal to zero. Returning to Equation 0, the estimate for m can be 
rewritten as 

m' = ((p^F-i)(FwF-i)(Fp))-i(p^F-i)(FwF-i)(Fcr), (6) 

using the fact that F^^F = 1. Finally, taking the Fourier transform of the various terms in 
Equation E] yields 

as an alternate expression to estimate the value of fh, where P = Fp, D = Fd, and W = 
FwF-^ Note that c = ( p-^wp) ^ = (P"^WP) ^ does not in general simplify as a result of 
Fourier transforming. 



5.2. The Bolocam Algorithm: Theory 

The science field maps produced by Bolocam each contain Hp ~ 20000 pixels, and an 
extremely large matrix must be inverted to calculate c since P^WP = p-^wp has dimensions 
of Hp X Hp. Direct inversion of such a matrix is possible, but is not practical on a typical 
high-end desktop computer. The map could be determined on a desktop computer via a 
conjugate gradient solver, but determining the covariance matrix would require a significant 
amount of simulation power. Therefore, we developed an algorithm to approximate rh' by 
exploiting the simplicity of our scan pattern, which involved raster scanning parallel to either 
the RA or dec axis. This approximation allows us to make maps in a relatively short amount 
of time using a standard desktop computer, which is extremely convenient. 

To illustrate this simplification, consider the map made from a single bolometer for 
a single scan within an observation. This scan will produce a one-dimensional map at a 
single dec value (for an RA scan) or a single RA value (for a dec scan). Each data point in 
the time-stream is separated by 24 arcseconds in map-space since our data are sampled at 
10 Hz and the telescope scans at 240 arcseconds/sec. Therefore, our data is approximately 
Nyquist sampled for Bolocam's ~ 60 arcsecond FWHM beams. The maps are binned with 
20 arcsecond pixels (1/3 of the beam FWHM, and slightly finer than Nyquist sampled), so 
p^ will map either one or zero time-stream samples to each map pixel. Note that Ug, the 
number of time-stream samples, will be slightly less than Up, the number of map-space pixels. 



^ Note that physical space values are denoted with a lower case letter, and the corresponding frequency 
space values are denoted with an upper case letter. 
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Since p"^ has dimensions of Us x rip, the sum of each row in is either one or zero and 
the sum of each column is one. Consequently, we will make the approximation that p^ = 1. 
From Equation El this means that 

(7) 

and therefore fh = d for a single scan of time-stream data. If we Fourier transform Equation[71 
then we find that 

C = w-\ 

Since W is diagonal, the inversion is trivial, and the result is that the Fourier transform of 
the map-space noise covariance matrix is diagonal with elements equal to the time-stream 
PSD*A/. 

The next step is to consider a map made from a single bolometer for a full observation, 
which contains twenty scans. We move the telescope in the orthogonal direction to the scan 
between scans by more than the size of a single map pixel, so we can still approximate p^ ^ 1. 
There are almost no correlations between scans because the atmospheric-noise subtraction 
coefficients are calculated scan-by-scan along with subtraction of the mean signal level. The 
covariance of maps made for a single observation from alternate scans is negligible, sup- 
porting this assumption that individual scans are uncorrelated. Therefore, the time-stream 
data and map-space data for different scans are essentially independent]^ Consequently, the 
noise in map-space will be stationary, which means that the noise covariance matrix can be 
diagonalized by Fourier transforming it. The Fourier transform of the full-map noise covari- 
ance matrix, C, can be visualized by noting that each diagonal element corresponds to a 
single Fourier-space map pixel (or equivalently, a single Fourier-space time-stream sample). 
So, this visualization of C will be equal to the single scan time-stream PSD* A/ for rows 
of map-space pixels that are parallel to the scan direction, and will have a white spectrum 
for columns of map-space pixels that are perpendicular to the scan direction. Alternatively, 
since there is a one-to-one correspondence between time-stream samples and map-space pix- 
els, this visualization of the diagonal elements of C is equal to the full map-space PSD*A/n, 
where A/q is the angular frequency resolution of the map. 

At this point, we need to add together all of the individual observations to make a single 
map. Since we have shown that the map-space data are equivalent to the time-stream data 
for a single observation, the easiest way to co-add data from separate observations is to use 



^ To verify that the data from different scans are independent, we created maps for each observation fr'om 
all the odd-numbered (right-going) scans and from all the even-numbered (left-going) scans. The cross PSDs 
of the right-going maps with the left-going maps were consistent with noise, indicating that the data from 
separate scans are independent. 
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the single observation maps. Since the noise in separate observations is uncorrelated, the 
maps can be co-added according to 

^=(^E^^^') Es-H' (8) 

where the subscripts i and j refer to observation number. The easiest way to evaluate 
Equation [8] is to Fourier transform it so that the noise covariance matrices are all diagonal. 
The result is 



M={j2^7'^ Y.^j'M„ (9) 

where M is the Fourier transform of the map and C is the Fourier transform of the noise 
covariance matrix, with diagonal elements equal to the PSD*A/n of the map. Since all of 
the Cs are diagonal, we can simplify Equation [9] to 

-1 



M 



where M is the two-dimensional Fourier transform of the map and V is the two-dimensional 
PSD of the noise in the map. At this point we have dropped the vector and matrix notation 
since Mj and Pi have the same dimensions. Note that A/q is the same for every map, so the 
constant factor of A/n from the first sum in Equation [10] cancels the factor of 1/Afn from 
the second sum in Equation [101 

Finally, to make a map using all of our data, we need to consider every bolometer, 
not just a single detector. To properly weight the data from each bolometer prior to co- 
adding, we calculate the expected variance, {<Jpf)j, in measuring the peak flux of a point-like 
source from a single scan through the center of the source for bolometer i. This variance 
is calculated using the scan-averaged time-stream PSD for each bolometer, PSDi{f), and 
the Fourier transform of the expected signal shape of a point-like astronomical signal, S{f), 
according to 

= ( f dfiMr^y' (11) 



PSD.if), 

where / is temporal frequency. Note that >S'(/) is the beam profile, not a delta function. 
Then, the data from each bolometer is weighted by a factor proportional to 1/(0"^^) prior to 
co-adding it with data from other bolometers. This is the optimal way to co-add the data 
for point-like signals; it is nearly optimal for signals of any shape if the PSDs have similar 
profiles for every bolometer, which is largely true for our bolometer signals since they are 
dominated by atmospheric noise. 



-17- 



However, due to atmospheric noise, along with our noise removal algorithms, there are 
correlations between the bolometers. But, most of these correlations are instantaneous in 
time and constant over the oteervation§ Additionally, the relative positions of the bolome- 
ters do not change during the observation, so the map-space separation of the correlations 
does not change. Therefore, the correlations are stationary in time with separations that 
are fixed in map-space, so the correlations are an additional time-independent covariance 
between map pixels that are sampled at the same time by different bolometers. This ad- 
ditional covariance is approximately stationary over the entire map, except where it breaks 
down near the edges because part of the focal plane is outside the map region. Since this ad- 
ditional covariance between map pixels is approximately stationary in space, its contribution 
to will be diagonal. Since is still diagonal, co-addition of the maps of individual 
observations can proceed according to Equation [TUJ Therefore, Equation [TU] can be used as 
the algorithm to produce our final science field maps. 

Note that we were forced to make several simplifying assumptions in order to develop 
Equation [TOl We have assumed that the pointing matrix, p-^, is equal to one. We have also 
assumed that the noise in our time-streams is stationary for each eight-minute-long obser- 
vation. Additionally, we have assumed that the PSD of the correlations between bolometers 
is white, and that all of the correlations are time-instantaneous. Finally, we have assumed 
that the map coverage {i.e., the number of time-stream samples that are binned in each 
map-space pixel) is uniform, so that the Fourier transform of the map is a valid description 
of the time-stream data. Deviations from these assumptions will alter the map estimate we 
compute from the optimal least squares map estimate. But, these deviations only affect how 
each time-stream sample is weighted before it is mapped. This means our final map will 
have more noise than an optimal map, but it will not be biased in any way. In other words, 
since the map-making operation is linear, the resulting map will be unbiased no matter what 
weightings are used to co-add the data, as long as the weights are properly normalized. We 
have confirmed this lack of bias via simulation, as we discuss below. 



5.3. The Bolocam Algorithm: Implementation 

To start, we must first produce a map from the time-stream data for each eight-minute- 
long observation. As mentioned in Section 15.21 this is done by calculating the variance in 



We have been able to find a small amount of correlated atmospheric signal that is not time-instantaneous. 
However, the time lag of these correlations is generally much less than one time sample, which means they 
will also be less than one map pixel. 
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measuring the peak flux of a point-like source under the assumption that the profile of the 
time-stream PSD is similar for every bolometer. To determine this variance, we calculate the 
PSD for each bolometer for each scan. These spectra are then averaged over all twenty scans 
for each bolometer, thereby making the assumption that the noise properties do not change 
over the course of the observation. Then, we determine the expected shape of a point-like 
source using our measured beam profile and scan speed. Finally, Equation [11] is used to 
determine the variance in measuring the peak flux of a point-like source for each bolometer, 
which is inversely proportional to the weighting factor applied to the time-stream data for 
that bolometer. 

At this point, we have individual observation maps for every observation, and we can 
make a map from all of the data using Equation [TDl But, one of the main assumptions made 
in developing Equation [10] was that the map coverage is uniform for each observation. If 
this assumption fails, then the Fourier transform of the map is not a good description of the 
time-stream data. Our scan strategy produced highly uniform coverage in the central region 
of the map, and this coverage falls rapidly to zero at the edges of the map. See Figure [6] 
To obtain sufficiently uniform coverage, we restrict our map to have sides of 42 arcminutes; 
the fractional RMS variations in coverage within this region for a single eight-minute-long 
observation are only about 8 - 9%. Since the coverage variations are minimal, we will assume 
that this square central region has uniform coverage, and therefore uniform noise properties. 
This assumption of uniform coverage allows us to directly compute the Fourier transform 
and noise properties of the map. We emphasize that, even if the assumption of uniform 
coverage fails and our algorithm is non-optimal, it is never biased because Equation [TD] is 
linear in the map-space maps. 

We now have a uniform coverage map for each observation, which can easily be Fourier 
transformed to produce the MjS needed in Equation [TU] But, we still need to determine the 
two-dimensional PSD of each single observation map. Due to residual correlations between 
bolometers, we do not understand the noise properties of our data well enough to determine 
the map PSD from simulation, so we instead estimate the PSD by generating a large num- 
ber of jackknifed maps from our real data. In each jackknifed map, a different subset of 
the time-streams from half of the scans within each observation was multiplied by —1. Note 
that the data from all of the bolometers within a single scan are multiplied by —1, so the 
residual atmospheric noise that is correlated between bolometers is preserved. This multi- 
plication leaves the noise properties of the map unchanged!"] while allowing us to produce 



Each time-stream sample (and therefore each map-space pixel) can be expressed as the sum of two 
signals: 1) an astronomical signal and 2) a random noise signal that is drawn from the underlying distribution 
of the noise in the Bolocam system. The astronomical signal corresponding to a particular map-space pixel 
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a large number of noise realizations for each map. Note that the residual atmospheric noise 
correlations are time-instantaneous, so they remain in the jackknifed realizations. We then 
generate 100 realizations for each observation, and we set the true PSD for each observation 
equal to the average of the map-space PSD computed for each realization. See Equation \T3[ 
Examples of the PSDs we calculated are given in Figure [3 This method of determining the 
map-space PSDs assumes that the time-stream data for each scan is uncorrelated with the 
data from all other scans, which we argued in Section 15. 2[ 

To determine the validity of the map-space PSDs we estimated from the jackknifed map 
realizations, we examined the distribution of PSD values for each realization. If the noise 
properties of the data are Gaussian, as we have assumed, then the PSD measured at any 
given Fourier map-space pixel will be drawn from 

f{X,,,) = (1/P,^)e(-^-^/^^^, (12) 

where Xj^j7 is the measured PSD for realization i at pixel u, Vjj is the true PSD for pixel 
V, and f{Xi^p) is the probability density function of Xi^p. Note that z7 has units of spatial 
frequency (i.e., radians"^), and describes a pixel in the spatial Fourier transform of the map. 
See Appendix |A] for a derivation of /(Xj ^;). The true PSD is estimated from 

i = Nr 

i=l 

where Nr = 100 is the number of realizations. To compare our measured PSDs to the 
probability density function (PDF) given in Equation [121 ^6 created the dimensionless value 

n.^ = fe, (14) 

with associated PDF 

f{Y,,) = e-^--. (15) 

Then, we compared our measured values of Yi^p to the PDF in Equation [151 In general, we 
found that our measured Yi^p follow a distribution extremely close to fiYi^p), except that 
the number of Yi^p with values near zero is slightly less than expected. Therefore, the map- 
space PSDs estimated from the jackknife realizations should be a good estimate of the true 
map-space PSDs. 



will be the same for any scan, and will disappear in the jackknife realizations when time-stream data from 
half of the scans is multiplied by —1. But, if the underlying distribution of the noise is Gaussian, then the 
distribution of signals it will produce is symmetric about 0. Therefore, the statistical properties of the noise 
will be unchanged when half of the data are multiplied by -1. 
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Let us consider the possible effects of imperfect signal removal in the single-observation 
jackknife maps. Because we only use the single-observation noise estimates as weights for 
coadding, the result of residual signal in the jackknife maps will be coaddition weights that 
are non-optimal. This non-optimality may degrade the noise of the final maps, but will not 
cause them to be biased. We may estimate the size of the signal leakage to determine how 
large the deviation from optimality could be. Our final fiat band power anisotropy upper 
limit is approximately 1000 /xK^^^. With our effective A ln(£) of 0.63 (derived in Section[6]), 
this upper limit corresponds to an excess variance in our final maps of ~ 500 /xK^^^^. These 
final maps have a variance of 10000 /iK^j^^^ (see Figure [9]). Given that ~ 500 observations 
contribute to each map, the single-observation map variance is ^ 5 x 10^ /xK^j^^j^, or ap- 
proximately 10000 times larger than our upper limit on the astronomical signal contribution. 
Even if we did not remove the astronomical signal using the jackknifing procedure, it would 
affect the single-observation PSDs, and therefore the weights, at only the 0.01% level. Using 
jackknife-generated PSDs reduces the effect of signal contamination further. Therefore, the 
effect of signal leakage into the single-observation jackknife maps is negligible. 



6. Transfer Functions 

The transfer function describes the fraction of the astronomical signal that remains af- 
ter processing as a function of map-space Fourier mode. In order to determine the transfer 
function of our data processing algorithms, we first generate a simulated map of the ex- 
pected astronomical signal. This map is then reverse-mapped into a time-stream using the 
pointing information in a real observation. Next, this simulated time-stream is added to the 
real bolometer time-streams from the observation, and then processed and mapped in the 
standard way. A map made from data that did not have a simulated signal added to it is 
then subtracted from this map, producing a map with the simulated signal after processing. 
Finally, the PSD of this map is divided by the PSD of the original simulated signal map to 
determine how much of the signal remains. Note that we are computing the transfer function 
for a PSD because we are interested in measuring an excess noise and not a specific signal 
shape, which means we do not need the phase of the transfer function. 

This transfer function was computed for twenty randomly selected observations, ten 
taken while scanning parallel to RA and ten taken while scanning parallel to dec. Realizations 
of the expected flat-band power anisotropy signal were used as the simulated signal. These 
realizations were generated in Fourier map-space assuming Gaussian fluctuations and a flat 
band power in Ce = Cei^d- + l)/27r of 50 fiK^j^^^. Note that although a fiat band power of 
50 fiK^j^j^ was used for the simulated signal maps, we found that the transfer function is 
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independent of the amplitude of the flat-band power anisotropy signal. For each observation, 
we averaged the transfer function obtained from 100 different signal realizations to determine 
the average transfer function. We then compared the average transfer function for each of the 
ten observations taken with a similar scan pattern. The result is that the transfer functions 
were the same within our measurement uncertainty for all of the observations. Therefore, 
we averaged the transfer function from all ten observations to produce a high signal-to-noise 
measurement for each atmospheric-noise removal method: average, planar, and quadratic. 
See Figure [HI 

Since all of the data processing is performed on the time-streams, the attenuation caused 
by the processing has a preferred orientation based on the scan strategy. The result is a trans- 
fer function that is not azimuthally symmetric because of the large amount of attenuation at 
low frequencies parallel to the scan direction due to atmospheric noise removal. Additionally, 
there is massive attenuation on scales larger than the Bolocam focal plane (~ 500 radians"^) 
because of the atmospheric noise removal algorithms. This occurs because these algorithms 
are designed to remove all time-instantaneous signals at each data sample, which is equivalent 
to subtracting any signals that vary slowly compared to the size of the focal plane. 

In addition to the signal attenuation caused by the data processing, the Bolocam system 
also attenuates some of the astronomical signal. By scanning across the sky, we are effectively 
convolving any signal with the profile of a Bolocam beam; since the beams have a non-zero 
width, this convolution will act like a low-pass filter on all of the astronomical signals. This 
filter will be approximately symmetric because the Bolocam beam profiles have a high degree 
of rotational symmetry. Additionally, since the beams are nearly Gaussian, the filter will 
be approximately Gaussian with a HWHM in variance of about 1000 radians"^ (which is 
equivalent to a HWHM^ ~ 6000 in angular multipole space). See Figure [HI 

In order to quantify the amount of signal attenuation by each atmospheric noise removal 
algorithm, it is useful to determine the effective bandwidth of the transfer function. The 
effective bandwidth describes the range of angular multipoles to which we are sensitive, 
as quantified by the transfer function, and can be used to convert an angular power, C^, 
to a map-space variance in /xK^^^^. In general, the effective bandwidth is calculated by 
integrating the transfer function over all angular multipoles. However, since the expected 
SZE power spectrum is approximately fiat in Ce, which results in a spectrum in Ce that falls 
like l/£{i + l), it is more useful to weight the transfer function by a factor of This 
weighting will produce an effective logarithmic, rather than linear, bandwidth, and can be 
used to convert an angular power in Ci to a map-space variance. This effective logarithmic 
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bandwidth, BWg//, is defined as 

BWe//= jduSlT^Bl 

J V 

where v is the two-dimensional spatial frequency, Si) is the expected signal spectrum, is 
the transfer function of the data processing in squared units, and B^} is the profile of the 
Bolocam beam. Since the expected anisotropy signal has a flat band power in Q, 



i{t + 1) 

for i = 27r|z7|@ Assuming this spectrum for S'i, a top-hat window between i = and 
i = irnax will producc a bandwidth approximately equal to 

BWe// oc \n{imax) - In(Cm) = Aln(^). 

Although the Bolocam transfer functions are not azimuthally symmetric, it is still useful to 
determine the effective A ln(£) for each of the atmospheric noise removal algorithms, with 
Aln(£) = 0.98, 0.58, and 0.37 for average, planar, and quadratic subtraction. Note that our 
final map, which consists of observations processed with different atmospheric noise removal 
algorithms as described in Section [3, has a bandwidth of Aln(£) ^ 0.63. 



7. Optimal Atmospheric Noise Subtraction 

Each of the science field observations were processed with average, planar, and quadratic 
sky subtraction, creating three separate files for each observation. Quadratic subtraction 
removes the most atmospheric noise, while average subtraction retains the most astronomical 
signal, so there is an optimal sky subtraction algorithm for each observation based on the 
type of astronomical signal we are looking for. To determine which algorithm is optimal, we 
computed a figure of merit, FOM, for each subtraction method. Since the anisotropy signal 
appears as a variance in the map, the variance on the amplitude of the anisotropy signal 
will be proportional to the square of the map PSD divided by the transfer function of the 
experiment. This can be seen in Equations IA7I and IA8I Therefore, the FOM is defined as the 
inverse of this variance on the anisotropy signal summed over all angular scales according to 

poM = ^i:?Il!g®!, (16) 



We have used the smah-scale flat sky approximation, £ = 2ti\v\. 
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where z7 is a two-dimensional spatial frequency with units of radians"^, S'l is the expected 
anisotropy power spectrum, is the transfer function of the data processing in squared 
units, -Bj7 is the profile of the Bolocam beam, and is the PSD of the noise in the map in 
squared units. Note that we have included the ~ 5 arcsecond uncertainty in our pointing 
model in Bp, and this pointing uncertainty effectively broadens the beam. To be precise. 



where Bp is the measured beam profile, and = l/27iap for a pointing uncertainty of 
For the anisotropy spectrum, we assumed a flat band power in Ce, so 



for i = 27r|z7|. The figure of merit is inversely proportional to the variance on an estimate 
of the anisotropy amplitude (in /iK^^j^), so it characterizes the signal-to-noise ratio of the 
map. 

In the end, average subtraction was the optimal method for just over 50% of the obser- 
vations, planar subtraction was the optimal method for just over 40% of the observations, 
and quadratic subtraction was the optimal method for just under 10% of the observations. 
We can calculate how much the observations optimally cleaned by each method contribute 
to our final S/N from 



where T denotes the set of observations optimally cleaned by a given method and FOMi is 
the figure of merit from Equation [161 The S/N contributed by the average/planar/quadratic 
observations is 70/29/1%. These ratios are different from the number of observations opti- 
mally cleaned by each method because the amount of atmospheric noise in the data generally 
determines which subtraction algorithm is optimal, and the observations optimally cleaned 
with average subtraction were made in the best observing conditions. Note that quadratic 
subtraction is the optimal method only when the weather conditions are extremely poor. 
This is because the anisotropy power spectrum falls quickly at high frequency, and the 
quadratic subtraction algorithm attenuates a large amount of signal at low frequency. For 
point-like sources, whose spectra are flatter, quadratic subtraction is the optimal processing 
method slightly more often. 



5, = B.e-I'^l'/^'^^ 




1 
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8. Final Map Properties 



Once the FOM is determined for eacli subtraction metliod for eacli observation, we can 
then produce a map of all of the data using the optimally processed map for each observation. 
To produce this final map, we need to make a slight modification to Equation [10] to account 
for the transfer function of the data processing and the Bolocam beam. We need to account 
for these effects because the transfer function depends on the scan direction and optimal sky 
subtraction algorithm for each observation. Therefore, the amount of astronomical signal in 
the map is in general different for each observation. To account for the amount of signal 
attenuation in each observation, the map PSD needs to be divided by the transfer function 
and the Fourier transform of the map needs to be divided by the square root of the transfer 
function. After making these modifications to Equation [TOl we have 



as the Fourier transform of the optimal map estimate, Ai. Tj is the transfer function of 
the data processing for observation i in squared units, Bi is the Bolocam beam profile for 
observation i, Mi is the Fourier transform of the map from observation i, and Vi is the 
noise PSD for observation i in squared units. Note that the astronomical signal in Ai will 
be equal to the true astronomical signal, because we have divided the Fourier transform of 
each single observation map, Mj, by the appropriate attenuation factor, a/tTB^H However, 
for some pixels in Fourier space, Tj and/or Bi take on extremely small values, which means 
that some pixels in both the numerator and denominator of Ai have extremely small values. 
Therefore, before taking the ratio of the numerator and denominator in Equation [T7] we 
apply a regularizing factor, so that 



M 




(17) 




(18) 



for 




(19) 



We have not included any phase information in the factor ^jTiBf because both the signal and the noise 
PSD contain only noise; the phase is irrelevant. 
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Although M' will be biased {i.e., it is not the Fourier transform of the true map of the sky), 
this bias is accounted for by the final transfer function we calculate in Section I8.2|[^ Note 
that M' can be Fourier transformed back to map-space to produce a map m', although m' 
will be biased. The maps, m', for each science field are given in Figure [9l 



8.1. Noise PSDs 

Analogous to the case of a single observation, we used jackknifed realizations of our data 
to estimate the noise PSD of m' . In this case, each realization is generated by multiplying 
a randomly selected set of half the observations in m' by —1. The map-space PSD from 
1000 realizations were averaged to determine the best estimate of the noise PSD for each 
science field, with the results shown in Figure [TOl In this section we establish that the noise 
PSD estimated in this way is statistically well-behaved (Gaussian) and unbiased. These 
characteristics are critical to the remainder of our analysis. 

We analyzed the distribution of individual realization PSDs to determine if the under- 
lying probability distribution describing the noise is Gaussian. As in the single observation 
case, we computed a dimensionless PSD value according to Equation [HI and compared the 
distribution of these values to the PDF given in Equation O In general, the agreement is 
good, indicating the underlying noise distribution is well approximated by a Gaussian. See 
Figure [T2l The Gaussianity of the noise PSDs of the jackknife maps is important because it 
justifies the form of the likelihood function we use. Equation IA5I (presented in Section [10. ip . 
That form assumes the Fourier coefficients of the final map are Gaussian-distributed random 
variables with variance given by the noise PSD estimated from the jackknife maps. 

Next, we show that the noise PSD estimated from the jackknifes is unbiased under 
two assumptions: 1) the covariance of any pair of distinct observations vanishes on average; 
and 2) negligible signal leaks into the jackknife maps. The first assumption is equivalent 
to the statement that there is no scan-synchronous or fixed-pattern noise in the maps. We 
have checked this assumption empirically and find that the average fractional covariance of 
distinct observations is ~ 2 x 10~^, which is consistent with noise. We will discuss below 



In Equations fTSl and fT9l TiBf /Vi acts as a weighting factor for each observation. Therefore, M 
represents the weighted mean of the Fourier transform of each single observation map divided by the square 
root of the transfer function for that map, (M/VtB^). Similarly, represents the weighted mean of the 
square root of the transfer function for each observation, VTB^. So, M' = {y/TB'^){{M/VTB'^)), which 
reduces to the weighted mean of all the single observation map Fourier transforms, M' ~ M, in the limit 
that all of the single observation transfer functions, TiBf, are the same. 
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how our non-detection of signal in the final map further justifies this assumption. With 
these assumptions, we can prove lack of bias of the noise estimate in a straightforward 
fashion by a simulation that obeys the assumptions. We generate a set of Nobs = 515 single- 
observation maps using the single-observation noise PSDs; these obviously have no signal and 
are uncorrelated with one another. Then, we construct Njack = 1000 simulated jackknife final 
maps and calculate the noise PSD of these maps, which also obviously have no signal. Next, 
we average these noise PSD estimates over all of the jackknife final maps, and divide by the 
input noise spectrum to determine how accurately we have recovered that input spectrum. 
The average (over all Fourier space pixels) of this normalized PSD is 0.9997±0.0006, showing 
that we indeed recover the input noise spectrum within the measurement uncertainty of our 
simulation. This exercise thus shows that the simulated final map noise PSD estimated by 
jackknife maps is an unbiased estimate of the simulated final map noise PSD. Note that we 
do not claim that the noise PSD generated in this fashion is the noise PSD of our true final 
map; the simulation is intended only to show that the jackknife noise PSD estimate method 
is unbiased. 

Let us now justify the assumption that negligible signal leaks into the jackknife final 
maps. In generating jackknife final maps, negative signs were applied to exactly one half of 
the observations. There are no fiuctuations allowed in the number of negative signs, only in 
which observations have them applied. Therefore, residual signal can arise in the jackknife 
final maps in only two ways: 1) if the relative calibrations of the different observations are 
imperfectly known or 2) if the weights of the different observations are unequal. 

In the first case, consider a single-observation fractional relative calibration error of ip, 
but assume all the component observations would otherwise be weighted equally (i.e., no 
variation in noise between observations). This fiuctuating relative calibration error does not 
cause a bias; in an ensemble of experiments, the final map has an expected signal value equal 
to the signal value of the true final map and the jackknife final map has an expected signal 
value of exactly zero. But, the calibration fiuctuations will cause an imperfect coaddition 
or cancellation of the signal in any given final map or jackknife final map realization, which 
will produce a fractional spread in the signal level of ip/y/ Nobs relative to the true signal. 
Given that ip < 3 /^c0 and Nobs — 500 observations per map, this error in both the final map 
and jackknife final maps is very small compared to the signal. Since we are not attempting 
a high signal-to-noise measurement, the vanishingly small size of the error relative to the 
signal is thus not a concern in the final map. The error affects the jackknifes in a more 



In Section [4?2l we calculated our flux calibration uncertainty to be approximately 5.5%. However, most of 
this uncertainty is due to systematics that will not change from one observation to the next. The uncertainty 
caused by fluctuations in the atmospheric opacity and the fit of our model are ~ 2 — 3%. 
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subtle way because it effectively adds noise to the jackl<:nife final maps, which means our 
noise estimate is slightly higher than the true noise level of our final map. However, this 
bias is negligible: the magnitude of the error is of order the signal times '^p/^JNohs■ We know 
ip < 3%, Nobs — 500, and the signal is less than 2% of the noise in the final maps (in RMS 
units, c.f. Section [53]) ■ so this bias is < 0.003% of the final map noise level in RMS units, or 
< 0.006% in variance units. The small size of the effect is not surprising: it is proportional 
to the signal size, and we have no detection of signal. Moreover, even if the effect were not 
negligible, it would result in an overestimate of the noise PSD and thus would result in an 
overly conservative upper limit. This kind of effect would only be problematic if signal were 
visible at high significance. 

In the second case, the argument is very similar, but now what matters is the fractional 
variation in observation weight. This fractional variation is large, < 2, due to the significant 
differences in atmospheric noise between observations. Here, the large number of observations 
and the fact that the spread is proportional to the signal size render the effect negligible. 
The RMS spread of the residual signal in the jackknifes will be < 2/ \/Nobs = 0.09 times the 
signal size. Again, because of the small size of the error relative to the signal and the lack 
of detected signal, the error due to this effect is insignificant. 

If the noise estimation approach has underestimated the noise (for example, by failing 
to account for non-stationarity or correlations, or in any other manner) then there will be 
more noise in the final map than expected from the jackknifes. However, because we find 
in Section 110.51 that our 90% CL interval on the amplitude of astronomical anisotropy. A, 
includes A = 0, we do not see any significant excess of noise in the true map above what 
is expected from the noise PSD estimate. This explicitly rules out scan-synchronous or 
fixed-pattern noise that would be averaged away in jackknife maps but would remain in the 
coadded map at the level of interest for this analysis. Had there been an excess above the 
expectation from the noise estimate, we would have had to show more explicitly that the 
noise estimate was correct in order to claim a detection. 

Alternatively, we consider the effect of overestimating the noise. The resulting final 
map PSD would be too low to be consistent with the noise PSD estimate. In our analysis 
(see Section [TU]) . this would yield a best fit value for the astronomical anisotropy A of zero. 
We do not find this to be true: since the best-fit value of A must lie inside the confidence 
interval of any value, and our 68% CL interval on A does not include A = (see Tabled]), 
the best-fit value of A must therefore differ from zero. 
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8.2. Astronomical Signal Attenuation 

Now that the noise properties of the maps are well described, we need to determine the 
amount of astronomical signal attenuation due to data processing, the Bolocam beam, and 
the regularizing factor in Equation [181 The method for calculating the transfer function of 
the data processing and regularizing factor is analogous to the method described in Section [H] 
for single observations. Contour plots of the total astronomical signal attenuation are given 
in Figure [TTl Compared to a single observation, the transfer functions for the final maps are 
much closer to being rotationally symmetric. The difference in the transfer functions is at low 
spatial frequencies parallel to either RA or dec, and is caused by adding observations made 
while scanning in perpendicular directions. This is because the modes in single observation 
maps, where there is a large amount of astronomical signal attenuation [i.e., at low frequency 
parallel to the scan direction), do not contribute much to the final map. Therefore, most 
of the signal at low frequency along the RA direction is obtained from maps made while 
scanning parallel to dec, and vice versa. This effect can be seen by comparing the plots in 
Figure El with the plots in Figure [TTJ 



9. Noise from Astronomical Sources 

Since the noise PSD of the final map is estimated from jackknifed realizations of the 
data, all of the astronomical signal will be absent from the noise PSD. This is fine for 
the anisotropy signal we are looking for, because we want to understand the noise of our 
system in the absence of our signal of interest. However, we need to estimate the amount of 
noise produced by sources other than the SZE-induced CMB anisotropies, including galactic 
dust emission, radio point-source emission, emission from dusty submillimeter galaxies, and 
primary CMB anisotropies. 

The amount of galactic dust emission can be estimated from maps of our science fields 



taken from the full-sky 100 fim DIRBE/IRAS dust map (IDIRBE websitd : ISchlegel et al. 



19981 ). To extrapolate the 100 /xm data to our band at 143 GHz ~ 2.1 mm, we have used the 



"model 8" extrapolation given in iFinkbeiner et al.l (119991 ). At 100 /xm, the typical surface 



brightness of the dust emission in our science fields is just over 1 MJy/ster, which corresponds 
to a surface brightness of around 5-15 bKcmb for Bolocam. Using the maps that have 
been converted to a thermodynamic temperature at 143 GHz, we determined the map-space 
PSD of the dust emission, which corresponds to a Ci less than 10"^ ^K^^^ for i > 10000 



Note that the resolution of the DIRBE/IRAS dust map is 6.1 arcminutes, which corresponds to HWHM 
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Since this is well below the expected SZE-induced CMB anisotropy we are looking for, it is 
safe to conclude that the signal from the dust emission in our maps is negligible. 

Emission from radio point sources will also contribute to the astronomical signal in our 
maps. The power spectrum from these sources can be calculated from 



Scut 



S^N{S)dS + wd'^ 



(20) 



where 5* is the flux of the source, N{S) is the differential number of sources at a given flux 
in a given solid angle. Scut is an estimate of the source-detection threshold in the map {i.e., 
the level at which sources may be detected and removed), Ct is the angular power spectrum, 
wi is the Legendre transform of the two-point correlation function of the sources, and 



SN{S)dS 



is the background contributed by the sources (iWhite and Maiumdarll2004l : iScott and White 
19991 ). We will assume wi = 0, since there is a large amount of uncertainty in the clustering 
of these sources]^ Di fferential number counts have been determ ined from measurements at 
1.4, 5, and 8.44 GHz JXoffolatti et al.lll998l : lOanese et aDll987h . with 



N{S) 



5GHz 



150 S-^-^ Jy-^ster-\ 



(21) 



Since the spectrum of the sources is nearly flat (i.e., oc z/^ with /9 = 0), this equation is 
valid over a wide range of frequencies. Additionally, the WMAP K, Ka, and Q bands have 



been used to determine the differential number counts at 22, 30, and 40 GHz (IBennett et al. 



20031 ). A^(>S') is similar for all three WMAP bands, and is < 70% of the value of the model 
in Equation [211 The differential number counts at 40 GHz are described by 



N{S) 



40GHz 



32 S-^'-' Jy-^ster 



To extrapol ate this equation to the Boloc am band center at 143 GHz, we will use the method 
described in lWhite and Majumdarl (120041 ). Since there is evidence of the power law for N{S) 



in £-space of < 2000. Therefore, we have no direct knowledge of the power spectrum on scales smaller than 
~ 6 arcminutes, which are the angular scales Bolocam is most sensitive to. However, the power spectrum of 
the dust falls rapidly at small angular scales, so the estimate a.t i < 2000 should provide a reasonable upper 
limit. 

Note that the total number of sources and total integrated power in i-space will not change if Wi is 
non-zero; the clustering modeled by wg will only shift power from high-^ to \ow-£. 
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flattening out at higher frequencies, they describe the differential number counts according 

iV(5)i43GHz = (20 - 32) 5^2.3 jy-igter-i. (22) 



We also need to estimate Scut in order to evaluate the power spectrum in Equation 
This cutoff flux will necessarily be somewhat arbitrary, but, since Ci is only weakly dependent 
on Scut, it will not significantly alter our result. We have chosen Scut = 10 mJy, which is 
approximately four times the RMS fluctuations per beam in maps made from data that have 
been optimally filtered for point sources^ Inserting this value of Scut into Equation [20} 
along with Equation [221 yields Q ~ 1.1 - 1.9 Jy^ steT~\ or ~ 7 - 12 x lO"*' fi'^luB- 
To compare this angular power spectrum to the expected SZE-induced CMB anisotropies, 
we determine the amplitude of a flat band power, C^^, that is required to cause the same 
temperature fluctuation as Ci given our transfer function, T^-Bf, according to 

^ 2i±lrp r)2 

^eff _ 47r ^^-"^ /9q\ 

For the radio point sources with = 7 — 12 x 10~^ /iK^jy^^, the effective Cn given the 
Bolocam transfer function is Cf^ ~ 35 — 60 /iK^^^, which is comparable to the expected 
signal from the SZE-induced CMB anisotropies. 

Additionally, emission from dusty submillimeter galaxies will be present in our maps. 
The same method used to determine the power spectrum from radio point sources can also 
be used to estimate the p ower spectrurn of these sources. We used the number counts 



distribution determined by lAguirrd (120081 ) , with 

iV(5)268GHz = 1619 5-2-26e-3035 Jy-^stcr^ 



The s pectrum of these objects can be described by 5*,^ oc , where 2.5 < /9 < 3.5 (IBorys et al. 



20031 ) . which gives a differential number count at 143 GHz of 

A^(5)l43GHz = (100 - 220) 5-2.26g-(2730-1460)S Jy-lg^er- 



There is some uncertainty in the spectrum of 5^ for these radio sources between 40 GHz and 143 GHz. 
White and Majumdar quote two spectra, one with /? = 0, and one with (5 = —0.3. This uncertainty in the 
spectrum of the radio point sources results in a finite range for the normahzation of the number counts after 
extrapolating to 143 GHz. 

From Equation [221 we only expect 1 — 2 sources brighter than 10 mJy in our entire survey of 1 
square degree, which is why we have not attempted to subtract out any sources prior to our anisotropy 
analysis. Additionally, the largest excursions in our maps are ~ 10 mJy, further justifying our choice to set 
Scut = 10 mJy. 
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Inserting the above formula into Equation [20] gives Ci = 0.4 — 1.2 Jy^ ster ^, or Ce 



3 — 9 X 10~® jiK.^(jjyiQ. Equation [231 can again be used to convert this to an effective constant 
Ci for our transfer function, giving Cf ~ 15 — 45 /xK^^jg. Ahernatively, we can compute a 
power spectrum usin g the differential number counts derived from SHADES data at 350 GHz 
(ICoppin et al.[ [20061 ). which is described by 



iV(^)350GHz = 2.2 X 10^ [S^ + (5.9 X lO^)^^-^]-^ Jy 



^^ster"^ 



Converting this N{S) to a differential number count at 143 GHz using the average spectrum 
of yields a similar power spectrum, with = 1.0 Jy^ ster~^, or (7/ = 8 x 10~^ uK 



which is consistent with the result from the number counts given by [Aguirrd (|2008[ ). 



2 

CMBi 



Finally, there will also be a signal in our map due to the primary CMB anisotropics, 
which are distinct from the SZE-induced anisotropics we are searching for. The power spec- 
trum of the primary CMB anisotropics has been measured to high precision by WM AP at 
i < 800 JNolta et al.ll2008h . and by ACBAR at 500 < ^ < 2500 jReichardt et al.ll2008h . This 
measured power spectrum is well fit by theory, with only a small number of free parameters. 
Therefore, we have generated a template of t he primary CMB power spectrum using the the- 



oretical prediction generated by CM BFAST (ISeljak and Zaldarriaga[[l996[ : [Zaldarriaga et al. 



1998[ : [Zaldarriaga and Seliaklboooh. with the best fit values to the free parameters from the 
WMAP 5-year data JPunklev et al.l[2008l ). Since the CMBFAST routine only computes the 
power spectrum up to £ = 3000, we fit a decaying exponential to the versus £ to ex- 
trapolate the primary CMB power spectrum to higher C.. We can again use Equation [23] 
to convert this power spectrum to an effective constant Ci given our transfer function, with 
Cf ^ 45 fiKl^.js. This band power is similar to what is expected from the SZE-induced 
CMB anisotropics. A summary of the expected signal from the various astronomical sources 
is given in Figure [T31 



10. Science Analysis 

10.1. Overview of Analyses 

In addition to instrumental noise from the bolometers, electronics, etc., our maps will 
contain an excess noise from astronomical sources, including anisotropics due to primary 
CMB fluctuations, fluctuations due to the SZE, and fluctuations due to unresolved astro- 
nomical point sources. It is our goal to constrain the level of these astronomically sourced 
noises, which we will specify as the amplitude of flat band power anisotropy power spectrum 
contributions in Ce. To obtain such a constraint, we must calculate the difference between 
the observed and expected power spectra of our maps and obtain a best estimate of the 
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excess noise, goodness-of-fit of the data to the model, including any possible excess noise, 
and confidence intervals for the amount of excess noise. This section describes how we obtain 
the estimate and intervals. 

The first analysis we perform will simply constrain the total astronomical anisotropy 
in the maps, without any interpretation of the source, assuming only that the astronomical 
noise has a spectral shape fiat in C^. 

The second analysis will statistically subtract the primary CMB anisotro py power spec- 



truna by using the precis e constraints placed on it by a variety of measurements (IReichardt et al 



20081 : iNolta et al.ll2008l ). The result will be a constraint on the non-primary-CMB contribu- 
tions to anisotropy, and will be mildly more sensitive because of the subtraction. We will do 
this by adding the expected "noise" from the primary CMB to our model of the instrumental 
noise. This expectation will fully take into account cosmic variance on the primary CMB 
anisotropy in a manner that we will explain below. 

In the end, this analysis will yield an upper limit on the astronomical noise. Because 
it yields an upper limit, it is conservative to immediately interpret the constraint as a limit 
on SZE anisotropy: if there are point source contributions, as we expect there are, then the 
SZE contribution will be smaller than the upper limit we obtain by the assumption that 
the point source contributions are negligible. The situation would of course be different, 
and that assumption would not be conservative, were we claiming a detection of excess 
non-primary-CMB anisotropy. 

One could extend this methodology to statistical subtraction of the non-negligible sub- 
millimeter and radio point source contributions, but the large uncertainties in those contri- 
butions as well as the possibly unknown systematic uncertainties lead us to conclude that 
the improvement in sensitivity will be negligible and somewhat untrustworthy. 



10.2. Deficiencies of a Bayesian Analysis 

We have chosen to model astronomical anisotropics using a fiat band power in C^, which 
corresponds to = AS'i, for £ = 27i\u\ and S'i = 27r/£(£ + 1). With these definitions, and 
assuming the noise PSD, Vp, fully describes the noise properties of the data for the reasons 
we have explained in Section [375| the best fit amplitude for an astronomical anisotropy signal 
is determined by maximizing Equation IA5t 



log{C) = V ( -\og{V, + ASlBlT,) - 



ueV 



V, + ASiBiT, 



with respect to A, where X^; is the measured PSD of the science field map in squared units, 
is the transfer function of our data processing in squared units, and is the profile of 
our beam@ Since our maps are real, X^t = X^jj, = V-.^, etc., so the sum only includes 
half the z7-space pixels, denoted by the set V . For reference, a detailed derivation of the 
above equation is given in Appendix |Al Note that Equation IA5I allows for A < 0. Although 
such values are not physical, fiuctuations in the noise can cause the most likely value of A to 
be less than zero when the expected value of A is small compared to the non-astronomical 
noise. 

The above expression is incorrect at some level because the i7-space pixels are slightly 
correlated, approximately 1 - 4% for nearest-neighbor pairs of pixels and less than 1% for all 
other pairs of pixels, while Equation IA5I treats all Fourier modes as independent. Note that 
the correlation function, c^^^i^ is largely translation invariant [c^^vi ~ civ — u')). This error 
due to pixel correlations raises three questions: 1) Does maximization of the likelihood given 
above result in an unbiased estimator of A7 2) Is this an approximately minimum variance 
estimator? and 3) Can we derive Bayesian credibility intervals on A from it? We have 
demonstrated using simulations that Equation IA5I remains an unbiased and approximately 
minimum variance estimator for A in spite of these correlations, presumably because ignoring 
these fairly uniform correlations does not shift the peak of C. See Table |2j However, the 
width of C is certainly dependent on these correlations: we are essentially over-counting the 
number of independent data points entering the likelihood and thus assuming more statistical 
power than we really have. 

We can make an approximate, unrigorous correction for the effective number of inde- 
pendent modes by calculating 



where A^true is the total number of z7-space pixels in V, Nes is the effective number of z7-space 
pixels, and Ci;y is the correlation between pixel u and pixel u'. The factor of 1/2 inside the 
parentheses arises from the fact that we have double counted the correlations with our sums 
over u and z7'; the factor of 1/Atrue is a normalization factor. c,7j7' is calculated from the 
Fourier transform of the map, M, according to 



Note that we are calculating the anisotropy amplitude for a single bin in ^-space. However, the technique 
can be applied to multiple bins in i-space by windowing the appropriate terms in Equation IA5I {i.e., if an 
£-space bin is described by the transfer function Tp, then IhVp, S'i T^vS];, and TpXp.) 




(|M,|)(|M,,|) 
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where the averages are taken over jackknife realizations. Equation IA5I is then multiphed by 
Ai'cfT / A^true to account for these correlations when calculating the Bayesian likelihood, with 
^cff/^true — 0.43 for our data. When ln(£) is exponentiated, this scaling factor will cause 
C to fall off less quickly than it would with A^^cfr = Atrue; thereby increasing the width of L. 

The lack of rigor behind the above correction implies that there will be statistical prob- 
lems in placing constraints using the above likelihood function. Were the above likelihood 
function correct, we could use it to set a a% Bayesian credibility interval on the parameter 
A by finding an interval [Ai, ^2], A^^A^ > 0, such that 

a J^yACjXM) 
100 J^dAC{Xp\A) 

where we have assumed a flat prior A > because of the non-physical nature of A < 0. If 
the likelihood function's width in A is not to be trusted, then such credibility intervals are 
not valid. Not even a simulation permits one to set a credibility interval because C is simply 
not the correct likelihood, even if its distribution can be characterized by simulation. 

Additionally, determining the goodness-of-fit for the best-fit value of A, A, will require 
simulation. That is, if Equation IA5I was correct, then we should be able to determine 
an analytic expression for the distribution of ln(£) for A that would allow us to calculate 
the goodness-of-fit of the data to the model. But, since the above likelihood function is 
incorrect, we must simulate an ensemble of measurements, with appropriate correlations in 
the Fourier modes, to determine the distribution of ln(£) for the value A. Therefore, the 
Bayesian approach offers no simplifications or reductions in computing time relative to the 
simulation-based frequentist technique we employ below. 



10.3. Overview of Frequentist, Feldman-Cousins Analysis Technique 

It is possible to deal with all of the above problems approximately correctly with a 
frequentist technique for establishing goodness-of-fit confidence levels and frequentist confi- 
dence (as opposed to Bayesian credibihty) intervals on A that incorporate the prescription s 



of Feldman and Cousins for dealing with a physical boundary (IFeldman and Cousins! 1 19981 ). 
The technique has two main features: 

1. First, we use jackknife maps with signal added based on an input value Agim in the 
physically allowed region Asim > to determine the distribution of C{Xi;^i{Asim)\Asim) 
as defined in Equation IA5I for an ensemble of experiments with outcomes Xj7 j(y4sim) 
With this distribution, we may determine whether C{Xff\Asim), the value of C for 
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the true data and the value Agim, is among the a% most hkely outcomes for that 
input value A^^^, thereby determining a goodness-of-fit confidence level. In doing this, 
we make the reasonable approximation that, although C is not a rigorously correct 
likelihood, it maps in a one-to-one, monotonic fashion to the true likelihood function 
£true- Specifically, if we consider two realizations Xi;^i and X,7^2, we assume that the sign 
of ln£(X,7^i|y4) — ln£(Xi;^2|^) is the same as that of In CtTue{Xu,i\A) — In £true(Xl;^2|^)• 
This assumption is far looser than the assumption that rescaling \nC by Xefj/iVtrue is 
correct; we are only assuming that the ordering of realizations in C and £true are the 
same, even if the numerical values are not the same. 

2. Second, we want to define a confidence interval of confidence level a% on A. Since 
we are taking a frequentist approach, these confidence intervals are defined to include 
the values of A for which, if A is the true value of the anisotropy amplitude, then 
the observed outcome is within the a% most likely outcomes (as defined below, a 
definition that is different than the usual likelihood) for that value of A. We use the 
same set of sim ulations with the follow i ng pro cedure based on the Neyman construction 



as modified by iFeldman and Cousins! (119981 ). We now calculate for each simulation 



realization i for each input parameter value Agim the ratio 

(X^^j (^sim) I ^sim) 



R^iA, 



^{^u,i{Asi^) \Ai 



where Asim is the simulation input parameter value (Agim > 0) and Ai is the best-fit 
value of A in the physically allowed region A > for the given realization i. We 
order the realizations in order of decreasing Ri{Asim) until a% of the realizations have 
been included; the value of Ri{Asim) defining this boundary is denoted by -Ra(Asim). 
The input parameter value Agim is then included in the a% confidence interval if the 
likelihood ratio for the real data, -Rdata(^sim), is among the a% largest Ri{Asim) values, 
-Rdata(^sim) > -Ra(^sim)- The interpretation is that, for values Agi^ belonging to the 
confidence interval of confidence level a%, the data is among the a% most likely 
outcomes, where "hkely" is quantified by -Rdata(^sim) instead of C{X0^i\Asira)- 

The above procedure can be conveniently visualized as follows. The simulations in- 
dicate that there is a smooth relationship between i?j(y4sim) and Ai at a given value 
of Asim- This is generically true, not specific to this analysis. Therefore, each simu- 
lation realization may be labeled by its value of A and we may write R{Asi^, A) in 
place of i?j(Asim)- We may visualize R{Asi^,A) as a function of A for a given value 
of Asim; the cutoff value Ra{Asij^) is a horizontal line in this plot, and so points with 
R{Asim,A) > Ra{Asira) map to a set of intervals in A; in fact, in our case, there is a 
single interval for each Asim- These intervals, called confidence belts, can be displayed 
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as intervals [y4i,y42] at the given value Agim in a plot of Agim vs. A, as illustrated in 
Figure [151 Do not confuse confidence belts, which are intervals along the A axis, with 
confidence intervals, which are intervals along the Asi^n axis as defined below. 

Then, to determine the confidence interval of confidence level a% on A given a data 
set Xi;, one finds Adata, draws a vertical line on the plot of confidence belts at the 
value v4data on the horizontal (A) axis, and includes all values of Agim for which the 
vertical line lies inside the confidence belt at that value of Asim- This confidence belt 
construction is equivalent to the above description based on R{Asi^, A) because the 
smooth relationship between R{Asim, A) and A ensures that, for a given Asim; if ^data 
is inside the confidence belt at a given value A^ira, then -Rdata > -Ra(^sim) for that value 

A . 



We comment on two important aspects of this construction of the confidence inter- 
vals. First is the ordering of the simulation realizations by R, not by £(Xj7_j(y4sim)|Asim)- 
Fel dman and Cousins discuss both possible constructions (the latter originally proposed 



by ICrow and Gardned (119591 )) and argue that the latter has a serious deficiency in that it 
ties the confidence level of the confidence interval to the goodness-of-fit confidence level; 
essentially, it is possible for the confidence interval to not provide the advertised frequentist 
coverage if the goodness-of-fit is poor. This typically happens when the experimental out- 
comes yield best-fit parameter values near or outside a physical boundary. In our application, 
this can occur if the simulation realization has a bit less anisotropy than expected, which 
would yield a best-fit A that is negative. In such a case, the (approximate) likelihood of the 
data set, C{X,;\Asitn), will in general be small. However, the approximate likelihood of that 
data set may not be small compared to the approximate likelihood, C{Xi;\0), of the most 
probable physically allowed alternative hypothesis of A, = 0. Feldman and Cousins show 
that, with this ordering principle, the confidence intervals never contain unphysical values 
for the observable. Additionally, there is a smooth transition from the case of an upper 
limit to a central confidence region, eliminating intervals that under-cover due to choosing 
between an upper limit and a central region based on the result. There is not room her e to 



reproduce their arguments in detail, we refer the reader to iFeldman and Cousins! (119981 ). 



The second important aspect is that the construction is done entirely by simulation so 
that the only way in which we depend on C, which we know to be deficient, is in the ordering 
it provides. We have assumed above that, in spite of its inaccuracy, C provides the same 
ordering of points as Arue, and hence A and R{Asi^, A) for a given (simulated or real) data 
realization and /^^(Asim) for a given Asim will be the same regardless of whether we use C or 

•^true ■ 
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10.4. Construction of Simulated Data Sets for Prequentist Technique 

To apply this method to our data, we first create a simulated map of the astronomical 
anisotropy for a given value of the astronomical anisotropy amplitude, Asim, using our as- 
sumed profile Si. This simulation is produced by drawing a value for each pixel, z7, from an 
underlying Gaussian distribution^ then multiplying it by AsimS^. The PSD of this simu- 
lated map is multiplied by our full transfer function and added to the jackknifed realization 
of our data, Xj ,70 Note that a different simulated map is created for each jackknifed real- 
ization of the data to allow for cosmic variance. Then, we use Equation I A5I to determine the 
most likely value of the astronomical anisotropy amplitude, Ai, for realization i. By using 
jackknifes of our actual data, we are including all of the correlations between pixels, and by 
simulating the astronomical anisotropy maps we are accounting for cosmic variance in the 
astronomical anisotropy. For a given value of Agim, we repeat this process for each jackknifed 
realization of the data. 

The data sets are then ordered based on the ratio of their likelihood to the likelihood of 
the most probable physically allowed outcome, i?j(Asim), as defined above. The procedure 
outlined in the previous section for defining Ra{Asim), finding confidence belts for each Asim, 
and then determining a confidence interval in Agim is then employed as described. 

To determine the goodness-of-fit of our data to the model given by 5^, we compare the 
likelihood of the actual data at the best fit value of Adata, 'C.{Xi;\Adata), to the likelihoods of a 
set of jackknifed realizations of our data with simulated spectra added according to with 
amplitude Adata, 'C{Xi;^i{Adata)\Adata)- FoT the observations of the Lynx field C^X^lAdata) 
is greater than C{Xp^i{Adata)\Adata) for 17% of the realizations, and for the observations of 
the SDSl field C{Xi;\Adata) is greater than C{Xi;^i{Adata)\Adata) for 43% of the realizations. 
Therefore, we can conclude that our model provides an adequate description of the data. 



We have also determined confidence intervals using non-Gaussian distributions for the SZE-induced 
anisotropy signal. The results are described in Section [10.61 

The reason we add the simulated astronomical anisotropy map to the jackknifed realization map instead 
of the time-streams is to reduce the amount of computational time required. Since the transfer functions 
of the maps are well measured, there is no reason to go all the way back to the time-streams to add the 
simulated signal. 
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10.5. Total Anisotropy Amplitude Results 

To determine tlie confidence intervals for the full data set, we make a joint estimate of 
A using both the Lynx and SDSl data sets. A plot of the Bayesian likelihood, along with 
confidence belts computed using the Feldman and Cousins method are given in Figure [T5l 
Uncertainties in our pointing model have already been included in these calculations by 
an effective broadening of the Bolocam beam. Our upper limits on the total anisotropy 
amplitude are equal to 590, 760, and 830 /iK^^j^ at confidence levels of 68%, 90%, and 95%. 
Note that the uncertainty on these limits due to the finite number of simulations we have 
run is ~ 10 — 15 /iK^j^^^. 

To determine the effective angular scale of our anisotropy amplitude measuremei its we 
have c omputed our band power window function, Wf using the method given by Knox 



( 119991 ). A plot of the peak- normalized band power window function for the full data set is 
given in Figure dni From this band power window function we have calculated an effective 
angular multipole for our data set, £eff, given by 

and equal to 5700. Additionally, the full- width half-maximum of the window function, 
FWHM^, is equal to 2800. A plot comparing our result to other measurements of the CMB 
on similar scales is shown in Figure [T71 



10.6. SZE-Induced CMB Anisotropy Results and Constraints on erg 

In order to determine the amplitude of the SZE-induced CMB power spectrum, we follow 
the same methods described above to determine the total amplitude of the anisotropy power 
spectrum. However, we now have to statistically subtract the signal due to the primary CMB 
anisotropics by accounting for both the amplitude and fluctuations of its expected power 
spectrum in the likelihood; these primary CMB anisotropics are effectively an additional 
noise in the map. The noise contributed to the map from the Bolocam system is given by 
Vp. Since the spectrum of the primary anisotropics in the CMB is well understood, we can 
calculate the expected noise from the primary CMB anisotropics. To calculate this noise we 
first create a simulated map of the primary CMB, assuming that the underlying distribution 



This band power window function is defined such that {Cb) — '^lO^f / ^)Ce, where (Cb) is the exper- 
imental band power measurement for the power spectrum, . Note that the transfer function of our data 
processing, Tg, is not the same as the band power window function, Wf /£. 
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of z7-space pixel values is Gaussian. This simulation is produced by drawing a value for 
each pixel, z7, from an underlying Gaussian distribution, then multiplying it by the best fit 
primary CMB spectrum given in Section O The PSD of this map is then multiplied by 
T^Bi and added to a jackknifed realization of our data, Xi^^, to give Xff^l A different 
simulated map is generated for each jackknifed realization of the data to account for the 
cosmic variance in the CMB spectrum. These modified jackknifed realizations of the data 
are then be used to determine the expected PSD, V^^^^\ for the noise contributed by the 
Bolocam system and the primary CMB anisotropics. We note that such a simulation of 
the primary CMB contribution is more correct than simply adding the primary CMB power 
spectrum to the non-astronomical noise power spectrum because it correctly reproduces 
pixelization and Fourier-mode correlation effects. 

Next, we select a model spectrum for the SZE anisotropics, S^^^^\ Using these new 
definitions, the Bayesian likelihood function in Equation IA5I can be written as 



[SZE] 



where A^^^^^ is the amplitude of the SZE-induced CMB anisotropics, Bi is the profile of 
our beam, and T^; is the transfer function of our data processing in squared units. As 
before, we create simulated SZE maps with an amplitude Af^^^\ add these to our jackknifed 
reahzations after mult i plying by {S]^^^^)'^B'iT^, then use the ordering method developed by 
Feldman and Cousins (11998 ) to determine the width of the confidence belt at A^f^^^\ By 
repeating this procedure for a range of physically allowed values of A^f^^^\ we can construct 
a full confidence belt that can be used to determine our confidence intervals. We emphasize 
that, while we use the Bayesian likelihood to construct our best estimator for A, a procedure 
that we have already demonstrated by simulation is unbiased and approximately minimum 
variance, we in no way rely on the Bayesian likelihood to determine confidence intervals on 
A. The frequentist Feldman-Cousins method is used for the latter task. 

Additionally, we need to account for the fiux calibration uncertainty. The uncertainty 
in the fiux calibration model derived from point sources is 5.5%, and the uncertainty in the 
area of our beam is 3.1%. Therefore, the uncertainty in our surface brightness calibration is 
6.3%. To determine the effect of this fiux calibration error on our confidence intervals, we 
multiplied each simulated primary and SZE-induced CMB map by 0j = 1 + where y is 
drawn from a Gaussian distribution with a standard deviation equal to our flux uncertainty 
of 0.063. A different (pi was generated for each simulated CMB map. This means that 
each simulated map has a different flux calibration, distributed according to our uncertainty 
in the calibration. New confidence belts were then calculated using the same procedure 
described above. We have also determined the confidence intervals assuming that there is 
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no uncertainty in the known flux of Uranus and Neptune [i.e., the only flux calibration 
uncertainties are due to our measurement errors and observational techniques). In this case, 
the flux calibration uncertainty is 3.5% instead of 6.3%. 

These flux calibration uncertainties produce non-negligible changes to the confidence 
intervals we determine for the anisotropy amplitude since it is a variance {i.e., it depends 
quadratically on the flux calibration). Therefore, for a simulated amplitude Agim, a fractional 
flux calibration uncertainty of cxj will increase/decrease the upper/lower bounds of the 68% 
CL confidence belt by an amount roughly proportional to Agimi^ + o"/)^- The resulting 
fractional change to the confidence interval limits will in general be non-trivial, but should be 
approximately equal to (I+ct/)^ for 68% CL limits. So, for a 3.5% flux calibration uncertainty 
we expect the 68/90/95% CL upper limits to increase by approximately 7/12/14% compared 
to the case of no flux calibration uncertainty. Similarly, for a 6.3% flux calibration uncertainty 
we expect the 68/90/95% CL upper limits to increase by approximately 13/22/27% compared 
to the case of no flux calibration uncertainty. After fully simulating the effect of the flux 
calibration uncertainty on our upper limits, we find results that are comparable to the 
predictions given above. See Table [H 

We have computed confidence intervals for two different SZE spectra: a flat spec- 
trum, (Sl-f^^^)"^ = 2 7r/i(i + 1) for i = 2tt\u\ and the analytic spectrum calculated by 



Komatsu and Seljakl (120021 ). The results for both of these spectra are similar, which is rea- 
sonable since the analytic spectrum is nearly flat at the scales to which we are most sensitive 
(4000 < ^ < 7000). See Tabled! In addition to the analytic spectrum calculated by Komatsu 
and Seljak, several SZE power spectra have been determined via hydrodynamic simulations 
using either MMH (moving-mesh hydrodynamic) or SPH (s moothed-par t icle hydrodynamic) 



algori t hms. Examples of MM H sim ulations can be found in IZhang et al.l (120021 ) , ISeljak et al. 



(l200l[ ). iRefrerier et al. 



2000f). and iRefregier and Tevssien (12002 ). Exa mples of SPH simu 



da Silva et al.l (120011 ) and ISpringel et al. 



(I2OOII ). Since most of the 



lations can be found in 
simulated SZE spectra are approximately flat at the angular scales we are most sensitive 
to, we have not determine d confidence levels using any of these spectra. See Figure 1 in 
Komatsu and Seljakl (120021 ). 



Komatsu and Seljak determined that the amplitude of the SZE-induced CMB anisotropies 
scale s according to aKflhhY an d is relatively insensitive to all other cosmological parame- 
ters (IKomatsu and Seljakl |2002| ). Using the results fror n the WMAP 5-year data, the best 
fit values for cxg, fit, and h are 0.796, 0.0440, and 0.719 (IDunkley et al.ll2008l ). These values 
produce a maximum SZE anisotropy amplitude of less than 10 fiK^j^^ at our band center of 
143 GHz for the analytic Komatsu and Seljak spectrum. For comparison, the 90% confidence 
level upper limit on the average value of the analytic spectrum weighted by the Bolocam 
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transfer function is 950 /iK^^g, including our flux calibration error. See Table [H Based 
on this upper limit, assuming the scaling relation given by Komatsu and Seljak and holding 
all other parameters fixed, the corresponding 90% confidence level upper limit on the three 
cosmological parameters is al{Qhhy < 2.13. Individually, the best constraint can be placed 
on as since the amplitude depends most strongly on this parameter, with < 1.55 at a 
confidence level of 90%. 

However, this upper limit has been derived by assuming the SZE-induced anisotropy 
signal is Gaussian, which is a poor assumption. To account for the non-Gaussianit y of the sig- 
nal, we have used a method similar to the one described by lGoldstein et al.l (I20031) to analyze 



White et al. 



Cooray 



data c ollec ted with ACBAR. Based on the results of numerical simulations by 
(2002) and lZhang et al.l (20021). a l ong w ith calculations of the trispectrum term from 
( 120011 ) and iKomatsu and SeljakI (120021 ) . they determined that the sample variance of the 
SZE-induced anisotropy signal should be a factor of three larger than the Gaussian equiva- 
lent for the £-range that ACBAR is most sensitive to. For our data, at £ ~ 6000, the sample 
variance is approximately four times larger than the expectation for a Gaussian. When we 
account for this increased sample variance our 68%, 90%, and 95% confidence level upper 
limits for the average amplitude of the Komatsu and Seljak spectrum are 790, 1060, and 
1080 /iK^jyj^, which are approximately 10% higher than the upper limits obtained from as- 
suming the fluctuations in the SZE anisotropy signal are Gaussian. The changes in the upper 
limits we determine are relatively minor because our uncertainty is dominated by Gaussian 
instrument noise rather than sample variance on the anisotropy signal. When we convert 
these upper limits on the anisotropy signal to an upper limit on ag, we find (Xg < 1.57 at a 
90% confidence level. 



11. Conclusions 

We have surveyed two science fields totaling one square degree with Bolocam at 2.1 mm 
to search for secondary CMB anisotropies caused by the Sunyaev-Zel'dovich effect. The fields 
are in the Lynx and Subaru/XMM SDSl fields. Our survey is sensitive to angular scales 
with an effective angular multipole of ieff = 5700 with FWHM^ = 2800 and has an angular 
resolution of 60 arcseconds FWHM. Our data provide no evidence for anisotropy. We are 
able to constrain the level of total astronomical anisotropy, modeled as a fiat band power 
in C£, with frequentist 68%, 90%, and 95% CL upper limits of 560, 760, and 830 fiK^jp^^. 
We statistically subtract the known contribution from primary CMB anisotropy, including 
cosmic variance, to obtain constraints on the SZE anisotropy contribution. Now including 
flux calibration uncertainty, our frequentist 68%, 90% and 95% CL upper hmits on a flat 
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band power in Ce are 690, 960, and 1000 (iKnMR - When we instead employ the analytic 



spectrum suggested by lKomatsu and SeljakI (120021 ) . and account for the non-Gaussianity of 



the SZE anisotropy signal, we obtain upper limits on the average amplitude of their spectrum 
weighted by our transfer function of 790, 1060, and 1080 fiK^^^. We obtain a 90% CL upper 
limit on as, which normalizes the power spectrum of density fluctuations, of 1.57. These 
are the first constraints on anisotropy and as from survey data at these angular scales at 
frequencies near 150 GHz. 

To calibrate the observations, beam maps were obtained using Uranus and Neptune. 
Pointing reconstruction was performed using frequent pointing observations of bright sources 
near our science fields. The data were flux-calibrated using technique s similar to those de - 



veloped to analyze earlier Bolocam survey data collected at 1.1 mm (iLaurent et al.ll2005l ). 
using Uranus and Neptune as absolute calibrators and a number of other sources as transfer 
calibrators. Internal uncertainty on the pointing and flux calibration contributes negligi- 
ble uncertainty to the final result; calibration uncertainty in the final result is dominated 
by uncertainty in models for the absolute brightness temperatures of Mars, Uranus, and 
Neptune. 

Our time-streams are dominated by fluctuations in atmospheric thermal emission (sky 
noise) and we developed several algorithms to subtract this noise from our data. We made use 
of our simple yet cross- linked scan strategy to develop a pseudo least-squares map-maker that 
can be run in moderate amounts of time on a single desktop computer. We used simulations 
to calibrate the transfer function of our data-taking and analysis pipeline and map-maker. 
We determined the expected noise properties of our final maps using jackknife realizations of 
the data obtained by randomly signed combinations of the ~500 independent observations 
contributing to each science field map. Our final confidence intervals on anisotropy level are 
determined using these jackknife realizations combined with the measured transfer function 
for anisotropies. 
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A. Appendix material 

The goal of our analysis is to determine the amplitude of the power spectrum due 
to emission from astronomical sources by measuring an excess noise in the maps of the 
science fields. This excess noise is the difference between the actual noise of the map, and 

the expected noise of the map based on measurements of the noise in the Bolocam system 
and knowledge of the expected signal spectrum. Therefore, we need measurements of the 
following quantities: 

• Xp: The measured PSD of the science field map at pixel P in units of ^K^j^^. V is 
a two-dimensional value, v — {I'RAj^dec), describing a location in the spatial Fourier 
transform of the map, and has units of 1/radians. 

• Vc^: The predicted PSD of the science field map at pixel u in the absence of the desired 
astronomical signal. Vt; is estimated from jackknife realizations, along with the PSDs 
of unwanted astronomical sources in the map {i.e., primary CMB anisotropies in our 
case) . 

• Sp-. The spatial power spectral profile of the expected astronomical signal. For a fiat 
band power (S| = 27r/ {£{£+!)), where the angular multipole i is described hy £ = 2n\iy\. 

• Bi: The peak-normalized square of the i^-space Bolocam beam profile. Since astro- 
nomical signals are attenuated by the beam, acts like a transfer function or filter. 
Note that the broadening of the beam in map-space due to our pointing uncertainty is 
included in B^. 

• Tp: The effective transfer function, or window function, of the data processing applied 
to the time-stream data. Analogous to Bi, Tp describes how much astronomical signal 
is attenuated. 
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With this convention, the expected PSD of the map can be described by 

{Xff)=V, + ASlBlT^, (Al) 

where A is the amphtude of the excess anisotropy power, in /iK^j^j^. 

The anisotropy amphtude can be estimated by determining what value of A maximizes 
the hkehhood of the measured map PSD, X^. Therefore, we need to determine the probabihty 
density function (PDF) describing Xj7, given A. First, note that 

Xi; =\a + ijSf , 

where a is the real part of the Fourier transform of the science field map and f3 is the 
imaginary part of the Fourier transform of the science field map. If we assume that the noise 
properties of the map are GaussianQ then the PDFs for a and /? are the same and are given 
by 

f(a) = ^^e-"'/'-' and f(p) = ^=e-^'/'''\ (A2) 

where o"^ = (Xp) /2. Next, after a change of variables to a = rcos(^) and (3 = rsin(6'), the 
PDF in Equation IA2I becomes 

Since the 6 dependence of /(r, 6) is trivial, we can reduce the above PDF to /(r) = 27r/(r, 9), 
with 

Finally, after one more change of variables using the relation Xc; = r^, we find that the PDF 
for Xj7 is equal to 

fix,) = :-^e-^'^/2'^^ (A3) 

where the factor of r has been replaced by 1/2 due to the change in the differential element. 
Equation IA3I can be written in terms of our measured parameters as 



This is an extremely good assumption. See Figure [T2l Although the anisotropy signal may not follow 
a Gaussian distribution, AS-BpTp <C Pp for a single v-sp&ce pixel, so the underlying distribution function 
for Xp will still be well approximated by a Gaussian. 
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using Equation IA1[ Note that we have made use of the fact that 2cr^ = (X,?) = Vci+AS'^B'iTij 
to go from Equation IA3I to Equation IA4[ 

The next step is to calculate a likelihood function, £, from Equation IA4I by multiplying 
f{X0\A) over all of the i7-space pixels. This product can be turned into a sum by taking the 
logarithm of £, with 



\og{C) = J2 [-^og{V, + ASjBjT,) 



V. + ASlBlT, 



(A5) 



Note that half of the i^-space is discarded from the sum in Equation IA5I since our maps are 
real {X^ = X_^). 

Then, the most probable value of the anisotropy amplitude for our measured map PSD 
can be determined by maximizing log(£) with respect to A. In practice, we maximize 
Equation IA5I by evaluating log(£) at a range of values for A. Since the number of i7-space 
pixels is < 10000, the computational time required to evaluate log(£) at each value of A is 
minimal, which means that we can determine the best fit value of A to almost any desired 
precision using this numerical method. 

However, it is also instructive to analytically approximate the value of A that maximizes 
Equation IA5I To start, we take the derivative of log(£) with respect to A, yielding 



d\og{C) 



dA 



= 0, (A6) 

A=A 



where 9,? = SiB'^T^ and A is the best fit value of A. For any given z7-space pixel, V,; ^ AQ^ 
for any physically reasonable value of A. Therefore, we can simplify Equation IA6I to 



~ 0. 

A=A 



If we rearrange some terms, and again keep only the lowest order terms in AQij/Vp, then we 
find 

A ^ V 



Finally, because AQp ^ Vjj, we can make the approximation that {X^) ~ Vp, which means 
that (2Xj; — Vv) — Vy. With this approximation we find 
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To understand this result, consider that for a single i7-space pixel the best estimate of A is 
{X^; — Vff)IQp. Therefore, Equation IA7I determines the weighted mean of A over all pixels, 
assuming that the uncertainty on the value of A for each z7-space pixel is proportional to 
'Pu/'^ui which is a reasonable assumption. This means that the variance on A implied by 
Equation IA7I is proportional to 

o\ oc (A8) 
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Table 1. SZE-induced CMB anisotropy results 



spectrum flux uncertainty 68% CL interval 90% CL interval 95% CL interval 



flat-total 





100 


- 590 iiKl^,j, 





- 760 iiKl^j, 





- 830 fxKlj,,^ 


flat-SZE 





90 


- 580 /iK^^5 





- 750 ^iKl^iB 





- 830 /iK^A.B 


flat-SZE 


3.5% (meas) 


90 


- 630 /iK2,^5 





- 790 fxKl,,^ 





- 880 fxKl,,^ 


flat-SZE 


6.3% (total) 


80 


- 690 /iK2,^5 





- 960 fiKl,,^ 


- 


- 1000 ^iKlj^s 


KS-SZE 





80 


- 540 ^lKl^B 





- 690 f^Kl^js 





- 770 ^iKl^^ 


KS-SZE 


3.5% (meas) 


80 


- 570 ^lKlMB 





- 740 fiKlj,,^ 





- 830 /iK^jvfB 


KS-SZE 


6.3% (total) 


70 


- 730 iiKl^^j^ 





- 950 fiKl^,^ 





- 990 fxKlj,,^ 


KS-SZE (nG) 


6.3% (total) 


90 


- 790 /iK^,,^ 


0- 


- 1060 /iK^,,5 


- 


- 1080 fiKl,,^ 



Note. — Confldence intervals for our estimates of the total and SZE-induced CMB anisotropy 
amplitude for both a flat SZE band power in Ce and the SZE spectrum given by the analytic 
model of Komatsu and Seljak (IKomatsu and Seljakll2002l ). The limits for the analytic model 
refer to the average amplitude of the SZE spectrum weighted by our transfer function. The three 
rows for each SZE spectrum give the upper limits for no uncertainty in our flux calibration, the 
3.5% uncertainty in our flux calibration due to measurement error, and the 6.3% uncertainty 
in our flux calibration due to the combination of measurement error and uncertainty in the 
surface brightness of Uranus and Neptune. The flnal row gives the confldence intervals when 
the non-Gaussianity of the SZE anisotropy signal is accounted for. 
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Table 2. Bias and efficiency of Equation IA5I 



input As 



average A 



real 



mm. var. a 



A 



I^^^MB 

100 ^iKl^B 
200 liKl^^j, 
400 ^iKl^jj, 
800 iiKImj, 



-2 /iK^jy.j^ 

96 iiKl^jB 
194 iiKl,,j, 
389 /.K^,,,^ 
806 /.K^,,,^ 



365 ^lKlJ^JB 

366 /iK^jv/B 

367 iiKl^jj^ 
371 /xK^,^^ 
380 /iK2,^,s 



12 jjK^j^^Q 
12 fiK^jj^^Q 
12 fiK(jj^^^ 
12 fiK(-,j^j^ 
12 fiK^j^Q 



270 /iK^,,,^ 
270 /iK2,^5 
270 /iK^,^^ 
270 /.K^,^^ 
270 /iK2,^,s 



Note. — A comparison between the amplitude of a simulated power 
spectrum added to a jackknifed realization of the data, Agim, and the 



most likely amplitude determined from Equation IA51 A. In each case 
1000 jackknifed realizations of the Lynx data were used, and the table 
lists the average value of A for these realizations along with the stan- 
dard deviation of the values of A. In each case the average value of 
A is consistent with A^im, indicating that Equation IA5I is an unbiased 
estimator of A. Additionally, to determine whether Equation IA5I is an 
efficient (minimum variance) estimator for A, we calculate the standard 
deviation of the estimates A for each input Agi^ and compare them to 
the standard deviation one would estimate using the Bayesian likeli- 
hood. Equation IA8I The latter underestimates the minimum possible 
standard deviation because the Bayesian likelihood is incorrect for the 
reasons presented in Section [l0.2[ Thus, the fact that the observed stan- 
dard deviation is only 40% larger than the Equation IA8l -based estimate 
gives us confidence that our estimator for A is reasonably efficient. 
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Fig. 1. — The solid black line represents a pre-down-sampled time-stream PSD, which has 
60 Hz pickup at frequencies above ~ 10 Hz. The dashed red line shows the time-stream PSD 
after downsampling and processing, including removal of an atmospheric noise template. 
The sharp increase in this PSD near 5 Hz is caused by the small amount of noise that is 
aliased in from frequencies just above 5 Hz during the downsampling process. Since there 
is approximately no astronomical signal at the frequencies where this noise increase occurs, 
this noise does not have a noticeable effect on our sensitivity. Overlaid as a dot-dashed green 
line is the beam profile, showing that very little astronomical signal will be present above a 
few Hz. 
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azimuth (arcsec) 



Fig. 2. — Location of the beam center of every detector relative to the center of the array. 
The red rings around each beam center represent the approximate FWHM of the beam. 
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Fig. 3. — All of the raw pointing data for Lynx at azimuth angles between -90 and 90 degrees. 
The pointing model (quadratic fit) is overlaid. Similar models were fit to the SDSl pointing 
data and the Lynx data at azimuth angles between 270 and 360 degrees. 
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8.1 8.3 8.5 8.7 

bolometer voltoge ot bios frequency (mV) 

Fig. 4. — Flux calibration for one of the six calibration data sets, overlaid with a linear fit 
of calibration versus bolometer operating resistance measured by the bolometer voltage at 
the bias frequency. Note that the bolometer voltage at the bias frequency is a monotonic 
function of the atmospheric opacity and the bolometer responsivity. 
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Fig. 5. — A histogram of the beam area calculated for each bolometer, with a Gaussian fit 
overlaid. 



single Lynx observation 



oil Lynx observations 




44°20m 



8h52m51m 50m 49m 48m 47m 
RA (J2000) 



o 
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o 




44°20m L 



8h52m51m 50m 49m 48m 47n 
RA (J2000) 



Fig. 6. — Map coverage, quantified by the number of time-stream samples tliat correspond 
to a particular map-space pixel for a single observation of the Lynx field and for the co- 
add of all observations of the Lynx field. The white square has sides of approximately 
42 arcminutes and contains the region of the map defined to have uniform coverage. The 
RMS deviations in coverage within this region relative to the average coverage within the 
region are approximately 8 - 9% for a single observation and around 1.5% for the co-add of 
all observations. 




Fig. 7. — Map-space PSDs, V^j, for single observations. The plot on the left shows an 
observation made in relatively good weather while scanning in the RA direction, and the 
plot on the right shows an observation made in relatively poor weather while scanning in the 
dec direction. In each case, note that there is a stripe of increased noise at low frequency 
along the scan direction, due to time-stream noise with a 1// spectrum. 
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average skysub quadratic skysub 




-2000-1000 1000 2000 -2000-1000 1000 2000 

RA wavenumber (rod"') RA wavenumber i/,^ (rod"') 



Fig. 8. — Contour plots of the transfer function, Tc^Bi, for observations made while scanning 
parallel to RA for average subtraction and quadratic subtraction. Note the large amount of 
attenuation at low frequencies along the scan direction and at scales larger than the focal 
plane size of approximately 500 radians"^. 
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Lynx SDS1 




ShSlm 50m 49m 48m 2h19m 18m 17m 

RA (J2000) RA (J2000) 



Fig. 9. — Maps of the science fields. Note that the astronomical signal in each map has been 
convolved with the transfer functions of the data processing and the beam, VtB^, but the 
noise has not been filtered in any way. The RMS of these unfiltered maps is approximately 
90 fJ^KcMB per beam, and the RMS after optimally filtering for point sources is ~ 70 fiKcMB 
per beam. 
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Lynx SDS1 




RA wavenumber (rod ') RA wavenumber (rod ') 



Fig. 10. — The map-space PSDs, Pci, of the maps made from co-adding all observations for a 
each science field. Note that relative to these PSDs, the power spectra of any astronomical 
signals will have been multiplied by the transfer functions of the data processing and the 
beam, T^B'i. 
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Fig. 11. — Contour plots showing the transfer functions, T^B'i, for the maps made from 
all observations of each science field. There is slightly more attenuation along the vra axis 
compared to the Vdec axis in the maps because more observations were taken while scanning 
parallel to RA compared to scanning parallel to dec. 
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Lynx 
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normalized PSD 



SDS1 




normolized PSD 



Fig. 12. — A comparison between the distribution of PSD values from the jackknifed reahza- 
tions to a Gaussian PDF for the data co-added over all observations for each Science field. 
See Equation [121 The agreement is good, indicating that the underlying noise distribution 
is approximately Gaussian. 
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angular multiplole (t) 



Fig. 13. — The power spectra from the primary CMB anisotropies (short green dashes), 
high and low estimates for radio point sources (red dash-dot), high and low estimates for 
submillimeter point sources (blu e dot-dot-dot-dash) , and t he analytically predicted SZE- 
indu ced CMB anisotropie s from iKomatsu and Seljald (120021 ) using the best fit value of as 
from iDawson et al.l (120061 ) (long orange dashes). Note that the point-source power spectra 
assume unclustered distributions. Also included as a solid black hne is the transfer function 
of the final map of the Lynx field with arbitrary normalization. 
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Lynx SDS1 




onisotropy amplitude (m^^^cmb) onisotropy amplitude (f^^cm) 

Fig. 14. — The Bayesian likelihood given by Equation IA5I for each science field. The likeli- 
hoods have all been normalized to one at the peak. These plots should only be considered 
as rough estimates for determining confidence intervals because the cosmic variance of the 
CMB spectra, correlations among map pixels, and the physical boundary that the anisotropy 
amplitude must be greater than or equal to zero have not been fully accounted for in the 
likelihood function. 
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Fig. 15. — The first plot shows the Bayesian hkehhood for a range of anisotropy amphtudes 
for the full data set, which includes all of the observations of both science fields. The 
remaining three plots show the frequentist confidence belts for the full data set for confidence 
levels of 68%, 90%, and 95%. 
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Fig. 16. — The band power window function, Wf /i, for the full data set. We have arbitrarily 
peak normalized the window function. 
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Fig. 17. — A plot of all of the current CMB anisotropy measurements above (. = 2000. Solid 
lines represent observations made near 150 GHz, and dashed lines represent observations 
made near 30 GHz. The primary CMB anisotropies are represented by a solid black line on 
the left side of the plot, and the predicted SZE-induced CMB anisotropies are shown as solid 
(150 GHz) and dashed (30 GHz) black lines. The primary CMB anisotropies were calculate d 
using the parameters given in Section [9l the analytic model of iKomatsu and SeljakI (120021 ) . 
along with the best estimate of from iDawson et al.l (120061 ) . were used to estimate the 
SZE-induced CMB anisotropies. All of the data are plotted with la error bars, except for 
the Bolocam upper limit at £ = 5700 and the BIMA upper limit at £ = 8748, which are given 
as 9 0% and 95% con f idence level upper limits, respectively. The ACBAR data w ere taken 
from iReichardt et al.l (120081). the BIMA data w ere taken from lDawson et al.l (120061 ). and the 
CBI data were taken from lMason et al.l (120031 ). 



